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Abstract 


The prediction of phase equilibrium behavior of complex mixtures over broad 
ranges of temperature and pressure is an important problem in chemical process 
design. Although cubic Equations Of State (EOS) are widely used in the prediction of 
Vapor-Liquid Equilibrium (VLE) data, they do not model the repulsive forces 
correctly and also they do not have the appropriate theoretical background. This work 
employs generalized Quartic EOS based on the Perturbed-hard-chain-theoiy, which is 
the basis of several successful EOS. The Quartic EOS incorporates dipole moment 
instead of critical pressure, which is used in cubic EOS. This work aims at prediction 
and comparison of VLE data for thirty-four different systems of seven different 
classes using four different models. Model- 1 uses the Peng-Robinson (PR) EOS with 
classical van der Waals one fluid mixing rule with single binary interaction parameter. 
Model-2 employs the PR EOS with composition dependent mixing rules of the 
Margules type with two binary interaction parameters, whereas Model-3 and Model-4 
relies on Qurtic EOS with the van der Waals one fluid mixing rule and composition 
dependent mixing rules, respectively. The results indicate that the Model-2 is the most 
suitable for the prediction of VLE data. It is also revealed that the mixing rules with 
two binary interaction parameters predict the VLE data much better than vdW single 
binary interaction parameter mixing rule. This is true for both PR and Quartic EOS. 
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CHAPTER 1 


EVTRODUCTION 

The estimation of phase equilibrium behavior of complex mixhires over broad ranges 
of temperature and pressure is an important problem in chemical process design as most 
of the chemical industrial processes deal with operations such as, extraction, adsorption, 
distillation, leaching, and absorption, which involve phase contacting. In particular 
distillation is the most important separation and purification process used in chemical 
industry, where vapor and liquid phases co-exist. This requires the process designer to 
have a reliable Vapor-Liquid Equilibria (VLE) data, for design of distillation columns. 
Though for a large number of binary and for a few multi-component systems the 
experimental VLE data can be found in literature, the experimental data is not available 
at all desired conditions. This necessitates the designer either to experimentally determine 
or estimate the required VLE data. But if the process of interest is to be performed over a 
wide range of pressures and temperatures, the number of experimental points required 
becomes large. Hence the available experimental data can be used to develop a model, 
which in turn can be used to predict the data at the required operating conditions. 

The available models for the prediction of VLE data can be broadly classified as 
activity coefficient models and Equation Of State (EOS) models. In activity coefficient 
model the liquid phase non-ideality is described in terms of an activity coefficient (yi) for 
component i of the system. The vapor phase non-ideality is expressed in terms of fugacity 
coefficient ((j)i) for component i of the system. The fugacity coefficient can be obtained by 
making use of an appropriate equation of state. Some of the models that belong to this 
method are Margules, van Laar, UNIFAC etc., (Rao Y. V. C. 1997). In EOS approach 
both the liquid and vapor phase non-idealities are explained by the fugacity coefficient 
where, it is indicated by (4^ for the liquid phase and (4h for the vapor phase. Some of 
the equations of state such as Peng-Robinson (PR), Soave-Redlich-Kwong (SRK), Peng- 
Robinson-Stryjek-Vera (PRSV), Benedict-Webb-Rubin (WBR) are used in this method 
(Mark et al. 1991). 



The activity coefficient models give satisfactory results at low to moderate pressures. 
However, their applications at critical and supercritical pressures yield inaccurate VLB 
data. EOS approach is a powerful tool for addressing this problem. The major advantage 
of the EOS approach is its applicability over a wide range of temperatures and pressures 
including critical and supercritical, where the activity coefficient models fail. In process 
simulation calculations, the equations of state are widely used due to their algebraic 
simplicity and their accuracy. 

There are two requirements for the equations of state to be successful. (1) The EOS 
must predict the saturation pressure of pure substances accurately. (2) The mixing rules 
must be available which correctly extend these equations to multicomponent mixtures. 
Different equations of state predict the thermodynamic properties with different degrees 
of accuracy and they require different types of data for proper evaluation. Hence, a 
reliable EOS for practical applications must be p>hysically soimd, reasonably simple and 
present no convergence problems. 

There has been significant improvement in the accuracy of cubic equations of state 
for the prediction of pure component properties in the last two decades (Sandler ct al., 
1994). Cubic equations of state, which are in wide use, model the repulsive forces 
^ incorrectly resulting in inaccuracies in representing the thermodynamic properties, 
especially in the supercritical fluid and condensed phase regions. Van der Waals 
proposed the first cubic equation of state in 1873. The well-documented van der Waals 
equation is the most influential EOS that incorporated both the repulsion and attraction of 
intermolecular forces in the expression. Numerous EOS have since been proposed within 
the firamework of this concept (Anderko (1990)). Redlich and Kwong (1949) modified 
the attractive term of the expression of the van der Waals EOS, but retained its repulsive 
term, to satisfy the boundary conditions in the low and high-density limits. The Redlich- 
Kwong (RK) EOS is among the most fimitfiil of the van der Waals family for routine 
engineering calculations. Various modifications of the RK equation to improve its 
accuracy have also been suggested. These are represented by the Soave (1972) and the 
Peng and Robinson (1976) equations. 



The modifications are mostly empirical and arbitrary, although in some instances 
physical constraints were proposed. All these equations of van der Waals family are 
explicit in pressure and third order in volume for which the equation can be solved 
analytically. The van der Waals EOS and its modifications account for the repulsive 
contribution to pressure with an excluded volume term RT/(V-b) that is conceptually 
incorrect at high densities. This term overestimates the repulsive pressure in comparison 
with the results from molecular simulations. The success of EOS of the van der Waal’s 
family is probably due to the mutual cancellation of errors in the repulsive and attractive 
terms. In order to overcome some of these deficiencies and to develop a theoretically 
based EOS which is capable of representing both the vapor and the liquid phases. Beret 
and Prausnitz (1975) proposed the perturbed-hard-chain-theory. This theory is the basis 
of several successful EOS. A feature of the perturbed-hard-chain-theory is the 
factorization of partition function into internal and external contributions. The density 
dependence of all external degrees of freedom is assumed to be same as those of the 
translational degrees of freedom, which is represented by Carnahan and Starling (1969) 
equation. Kubic (1986) proposed the Quartic hard chain equation of state. Shah et al., 
(1994) proposed the generalized Quartic EOS for pure non-polar fluids. Later Shah et al., 
(1996) proposed the Extended Generalized Quartic equation of state for pure polar fluids. 
The Quartic EOS retains the desirable features of the popular cubic EOS with improved 
accuracy. 

The extension of equations of state to describe phase behavior for a broad range of 
multicomponent mixtures had been more difficult due to limited applicability of the van 
der Waals (vdW) one-fluid mixing rules that are commonly used for relatively simple 
mixtures. The vdW mixing rules are incapable of representing the highly non-ideal 
mixture behavior of polar or associating fluids. Consequently, much effort in recent years 
has been devoted toward developing alternative mixing rules to extend the cubic 
equations of state to complex mixtures [Orbey and sandier (1996); Heidemann (1996)]. 
With the advent of multiparameter mixing rules, equations of state are now being used 
for the phase equilibrium calculations of complex mixtures that were traditionally 
described with activity coefficient models. 



The present work deals with a comparison of Quartic EOS with Peng-Robinson 
(PR) EOS for the estimation of VLE data using two different mixing rules. The 
generalized quartic equation of state has been used to predict the VLE data. The quartic 
EOS has four roots. The identification of these roots is very easy as, one root is always 
negative which has no physical significance. Hence, the rules used to identify the roots 
of cubic EOS could be used with quartic EOS also. This makes the prediction of volumes 
at different temperatures and pressures possible, by solving the corresponding equation of 
state in terms of the compressibility factor (Z) or volume (v), 

VLE calculations are done using Quartic EOS and PR EOS, for binary nuxtures 
of hydrocarbons, alcohols, and refiigerants. Van der Waals one fluid mixing rules and 
composition dependent mixing rules [Huron and Vidal (1979); Stryjek and Vera (1986)] 
are used with each EOS for the VLE calculations. The results are then compared with the 
available experimental data. The Quartic EOS with composition dependent mixing rules 
predicts accurate VLE data for alcohol-alcohol mixtures. For alcohol-hydrocarbon 
mixtures the PR EOS with composition dependent mixing rules predicts accurate VLE 
data except for the mixtures that contain butanol as one component for which the Quartic 
EOS with composition dependent mixing rales predicts accurate VLE data. For 
hydrocarbon-hydrocarbon mixtures the PR EOS with composition dependent mixing 
' rales is as good as the Quartic EOS with composition dependent mixing rales. For 
reJfrigerant mixtures, in most cases, the PR EOS with composition dependent mixing rales 
predicts extremely accurate VLE data. However, the Quartic EOS with composition 
dependent mixing rales is also equally good. 



CHAPTER 2 


PHASE EQUILIBRIUM AND MIXING RULES 

2,1 Equilibrium Criterion 

For a closed system containing vapor and liquid phases, with each phase consisting 
of c components, in a state of equilihrium at constant temperature (J) and pressure (P), 
the criterion for equilibrium is given by (Smith et al. 1996), 

= c) (2.1) 

where, /r/ and nj are chemical potentials of component i in liquid and vapor phases, 

respectively. The chemical potential ju,. , is given 

Ai,=r,(j’)+i?nn/;. ( 2 . 2 ) 

A 

where T,. (T) , an integration constant, is a function of temperature only and is the 

fugacity of component i in the mixture. Using Eq. (2.2) in Eq. (2.1), the general criterion 
for vapor-liquid equilibrium reduces to 

/•=// (^- = 1.2, c) (2.3) 

The Eq. (2.3) can be rewritten as (Rao Y. V. C. 1997), 

0=1.2 c) (2.4) 

where y,- = activity coefficient of component i in the solution. 

Xj, yi = mole fractions of component i in liquid and vapor phases, respectively. 

= fugacity of component i in the standard state. 

A 

(/>/ = fugacity coefficient of component i in the vapor phase. 

P = Pressure at which the system is held. 

The standard state fugacity {f°) is the fugacity of component i at the mixture 

temperature and at some arbitrarily chosen pressure and composition. The standard atate 
for component i, is chosen as the liquid phase of the pure component i at the system 
temperature and pressure. With this choice of the standard state, is given by. 



= P/0/ exp 


(2.5) 


f v/(p - p^y 

1 j 

where, is the fiigacity coefficient of component / at saturation pressure f*‘ for the 

same component i at temperature T. v\ is the molar volume of liquid for component i. 

Substitution of Eq. (2.5) in Eq. (2.4) results in the basic equation for vapor-liquid 
equilibrium, which is given by 

(; = i,2, c) (2.6) 

At low pressures, the vapor phase can be assumed to behave like an ideal gas. Hence 

A 

(f)!' =1 and </>/=!. Moreover, the Pointing correction factor (exponential term in Eq. (2.6)) 

A 

is approximately equal to unity. At moderate pressures and (^/ are approximately 


equal and hence it is reasonable to assume that (<f>J /((>■) - 1- Thus at low to moderate 
pressures Eq. (2.6) reduces to 


YiXiP- =yiP 

(1,2 


(2.7) 

II 

(1,2 

c) 

(2.8) 


' However, at high pressures these assumptions are not valid. Eq (2.3) can also be written 
as 


or 


or 


i‘x,P = ijy,P 

A 

(since = 

(hhiP ) ) 


ilx,=ijy, 

(1,2 

c) 

(2.9) 

II 

=^1 

II 

(1,2 

■c) 

(2.10) 




where is the K factor or equilibrium constant for component i. 

The methods available for the prediction of VLE data can be broadly classified as 

1. Activity coefficient Approach 

2. Equation of State (EOS) Approach 



The main difference between the two methods lies in the reference state. For the 
EOS method, the reference state is the ideal gas condition and the fugacity coefficient is a 
measure of the deviation from this state for both the vapor and liquid phases. For activity 
coefficient methods, the reference state for the liquid is the pure liquid at the same 
temperature and pressure of the mixture. The choice of the thermodynamic model to be 
used depends on 

• Nature of the components 

• The pressure and temperature conditions 

NATURE OF THE COMPONENTS 

The nature of the component decides the extent of non-ideality. Non-ideality is due 
to the interaction between the molecules, which depends on intermolecular forces. This 
interaction is more evident in the liquid phase where the molecules are packed closer. 
The choice between the EOS method and Activity coefficient method can be made based 
on the polar nature of the components. 

Hydrocarbons like ethane, propane, pentane are non-polar and can be represented 
adequately by EOS methods. Mixtures of polar compounds may be represented using 
Activity coefficient methods. Mixtures of polar and non-polar compounds can give rise 
' to two immiscible liquid phases and hence a model that is capable of predicting a phase 
split must be considered. This would be an Activity Coefficient method like NRTL or 
UNIQUAC. Wilson or Van Laar models are incapable of predicting splitting of the liquid 
phase. 

TEMPERA TUE AND PRESSURE 

At high pressures and at temperatures above critical, the activity coefficient models 
fail. EOS models are valid over a wide range of temperature and pressure. Since 
consistency is inherent in the EOS approach, these are applicable even at the critical 
conditions. 



2,2 Activity Coefficient Approach 

The basis of activity coefficient approach is the activity, which is defined as 

a - (2 10 

/,nr.rx,-) 

Here a^ is the activity of component i and f° is the standard state fugacity (at some 
arbitrarily defined pressure P° and some arbitrarily defined composition x°). The 
activity coefficient y,. is defined as 


y. =-^ = . 


The partial molar Gibbs free energy is given by 


( 2 . 12 ) 


E 

M ° = gi.real “ S i, ideal = = RT\Xi- ‘ 


fi,i. 


.ideal 


From Roults law one can obtain 


fmeal =Xif° 


tO 

g‘ =RT\RlJ^^ = RT\Rr, 

Using the definition of partial molar excess Gibbs free energy, one can obtain 


(2.13) 

(2.14) 

(2.15) 


Iny,. 


RT 






(2.16) 


' ^T,P,nj 


fr • p m m 

where, g = molar excess Gibbs energy. The relation between g and the activity 
coefficients can also be expressed as 


= T^igi^ = RT'Zxi Iny,. 


/=1 




(2.17) 


The Gibbs-Duhem relation for a binary mixture is given by 


5 Iny, 


V 5^1 JrP 


^51ny. ^ 

V ^-^2 Jt.p 


(2.18) 


Assume that f can be expressed as a function of temperature, pressure and composition. 
That is. 


g^=f(T,P,x^) 


(2.19) 



The Activity coefficient can then be derived from Eq. (2.16) and substituted into the 
following equilibrium constraint. 

yiP^J =^yi(T,p,x,)xj;‘ ( 2 . 20 ) 

The system of equations defined by Eq. (2.20) can be solved for the equilibrium 
temperature, pressure, and concentration. Therefore, the basis of activity coefficient 
approach is to find an accurate form of Eq. (2.19). 

The commonly used Activity coefficient models are 

> Margules Equation 

> Van Laar Equation 

> Wilson Equation 

> NRTL Equation 

> UNIQUAC Equation 

> UMFAC Equation 


2.2.1 Margules Equation 


The Margules equations, can be deduced from the expressions for excess Gibbs free 
energy given by WohTs expansion (Mark et al. 1991) 


(< 




RT{x,q^+x^q^) 


— 2Uj 2.Z,Z2 + 3ai],2X, Z 2 + 30,22-2122 + 


where 


( 2 . 21 ) 


X.Q . 

2 . = ili = the volume fraction of component i ( 2 . 22 ) 

x,q, +X2^2 

qi = measure of volume of the component i. a^, aij 2 , aj 22 are empirical constants 
which signify two body, three body,. . .. interactions. 


Margules assumed the size of molecules to be equal {qi = q 2 ). After some re- 
arrangements Eq. (2.21) reduces to the Margules three suffix equation given by 

~ (-^12 -^21 ^2 )^ 1'^2 
K1 

where An and ,^ 7 ; are adjustable parameters. 


(2.23) 



The activity coefficients y, and are given by 

^Y,=lA2+2iA„-A„)x,}cl (2.24) 

Iny 2 = ~ -^21 )-^2 ]^i (2.25) 

The Margules three suffix equations are used for symmetrical systems (qj = qi) i.e. 
systems where all the components have almost equal size. 


2.2.2 Van Laar Equation 


Margules equation assumes equal size for molecules, which is not true. If the size 
difference of the molecules (^, ^qj) is incorporated into WohTs expansion, (Mark et al. 
1991) it gives 


RT x,^, + ^2^2 

Simplifying Eq. (2.26) leads to van Laar equations, which are given by 

g A^2-^2lXiX2 


RT x^Aj2+X2A2^ 


Iny, = A , 2 


lnY2 = A2 



_XjAj2 +X2A21 

where A n and A 21 are adjustable parameters. 


(2.26) 


(2.27) 

(2.28) 

(2.29) 


2.2.3 Wilson Equation 

For Mixtures in which the components differ from each other in molecular size and 
interaction between the xmlike and like molecules are different, Wilson (1964) proposed 
the following expression for the excess Gibbs free energy of a binary solution. 

E 

^ = -X, ln(x, + A 12 X 2 )- X 2 ln(x 2 + Ajix, ) (2.30) 

where An and A 21 are two adjustable parameters, which are related to pure component 
molar volumes and characteristic energy differences. 



Ai 2 and A21 are given by 


A ^2 

Aj 2 = — exp 


Aj, =-!-exp 


^12 ~^21 

RT 


A-J2 A,2i 


RT 


V, —molar volume of pure component i. 

A-ij=energy of interaction between molecules of component i and j. 

The activity coefficients can be derived from Equation (2.30) and are given by 

Iny, =-ln(x, +A, 2 Ji: 2 )+x 
lny 2 = ~ln(x 2 + A21X, )— Xj 

X, + AjjXj Xj + AjjX, 



(2.31) 

(2.32) 


(2.33) 

(2.34) 


2.2.4 Non-Random Two Liquid (NRTL) Equation 


The Non-Random Two Liquid (NRTL) equation proposed by Renon and prausnitz 
(1968) is applicable to partially miscible as well as completely miscible systems. 
The NRTL equation for excess Gibbs free energy is given by 

(2.35) 

where 


g 


RT 


:X,X 2 


^21 ^21 ^ 12^12 


X, + 0,2X2 X2 + G 2 ,X, J 


_ gu gzi 
RT 


(2.36) 


_ Six 

21 

(2.37) 

<^12 = exp(-a,2T,2) 

(2.38) 

O 2 , =exp(-a2,T2,) 

(2.39) 


g-y—energy of interaction between molecules of component i and j. 

The parameter an is introduced to take into account the non-randomness of the mixture. 
If an=0, the mixture is completely random and Equation (2.35) reduces to Margules 
equation. The activity coefficients can be derived from Equation (2.35) and are given by 



l^y ^ ^21 . M2^12 

I Xj + ^ 21^7 1 (^2 ■*■ ^ 21^1 ) 


■^12^12 


In/j = X, T 


’^ 21^12 


X2+Gj2xJ (jCi+Gi2X2)^ 


(2.40) 


(2.41) 


The NRTL equation does not provide any additional advantage over Margules and Van 
Laar equations for moderately non-ideal solutions, but it represents the excess Gibbs free 
energy of strong non-ideal and partially miscible systems quite satisfactorily. 

2.2.5 UNIversal QUAsi Chemical (UNIQUAC) Equation 


To express the excess Gibbs free energy of a binary mixture, Abrams and Prausnitz 
(1975) developed the UNIversal QUAsi Chemical (UNIQUAC) model. The UNIQUAC 
equation for (^/RT) contains two parts - a combinatorial part and a residual part. The 
combinatorial part takes into account the composition, size and shape of the constituent 
molecules and contains pure component properties only. The residual part takes into 
account the intermolecular forces and contains two adjustable parameters. The 
UNIQUAC equation is given by 


where 


— — = — — (comhinatoridS) ■¥ — — (residual) 
RT RT RT 


A . z , ^1 1 ^2 

-= — (combinatorial) = x. In— ^ + X 2 In 1- — In h ^ 2^2 

RT X, X2 2 i ^2 


JC> 

^(residual) = -q^x^ ln(0, + 02 ^ 12 ) “^ 2^2 +^ 1 ^ 21 ) 

RT 


(^/=segment fraction of component i = Y.Xjrj 


0/=area fraction of component z= x^qJj^Xjqj 


(2.42) 


(2.43) 


(2.44) 


(2.45) 


(2.46) 


r,- = volume parameter of component i 
qi = surface area parameter of component / 



T.. = exp' 


=adjustable parameter 


(2.47) 




RT 


exp 


V 


Uij=svQX3%G interaction energy for the interaction of molecules of component i with the 
molecules of component and 
z = coordination number which is usually taken as 70. 

From Eqs. 2.42-2.47 one can obtain the activity coefficients as 


In/,. = \xxY ipombinatorial) + lny,^(ras/dwa/) 
Iny/ =ln^ + ^^,. ln^ + /, 

X, 2 0,. 
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( 2 . 49 ) 

(2.50) 


where /. = -q,)- {r. -l) (2.51) 

The structural parameters n and q,- are calculated as the sum of product of the group 
volume and parameters 72* and Qk. 

(2.52) 


?,=Ev''^a (2-53) 

k 

where v/'^ is the number of groups of type A: in a molecule of component i. The 
UNIQUAC equation is applicable to a wide variety of liquid solutions commonly 
encountered by chemical engineers. 


2.2.6 UNIquac Functional group Activity Coefficients (UNIFAC) Method 
(group contribution method) 

To correlate and to predict the thermodynamic properties of a solution it was 
convenient to consider the component molecules as a collection of functional groups. 
The functional groups are structural units such as -CH 3 , -OH and others which when 
added form the constituent molecules. In the group contribution methods, a solution of 
components is treated as a solution of groups. The activity coefficients of the 


components are then determined by the properties of the functional groups rather than by 
those of the molecules. The basic idea behind this approach is that there are several 
thousands of chemical compounds that are of interest to a chemical engineer but the 
number of functional groups, which constitute these compounds, is very small. The 
activity coefficients in a large number of solutions can be calculated from the parameters 
characteristic of the functional groups. 


In the UNIFAC method the activity coefficients consist of two parts, a combinatorial 
part ) and a residual part {Inyt ). The combinatorial part of activity coefficient Inyf^ 
given by Eq. (2.46) of UNIQUAC method is directly used. To evaluate Inyf only pure 
component data is required. One identifies the functional subgroups present in each 
molecule and estimates the parameters r,- and qi as the sum of the group volume and area 
parameter Rk and Qk by using the Eqs. (2.52) and (2.53). 

The residual part of the activity coefficient /«y/ is given by 

toy," =pl'’(tor, -tory>) (2.54) 

where 
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(2.56) 


RT j 

Xm = mole fraction of group m in the mixture 

amn = group interaction parameter (in K). ^mn ^ ^nm 

Um„ = measure of interaction energy between groups m and « 

Fk = group residual activity coefficient; and 

= residual activity coefficient of group A: in a reference solution containing 
only molecules of type i. 


(2.57) 


2.2 . 7 Advantages and Disadvantages of Activity Coefficient Approach 

A simplifying assumption that is made in the development of activity coefficient 
models is that is independent of pressure. In view of this assumption, the activity 
coefficient approach can be applied to systems at low to moderate pressures, far removed 
from the critical. For processes that operate near or above the critical, the activity 
coefficient approach cannot be used, which is by far the greatest deficiency of this 
approach. 

The activity coefficient approach should be used to treat low pressure VX,E data. The 
activity coefficient approach can be successfully applied to various types of systems: 
mixtures of alcohols with hydrocarbons, non-polar hydrocarbon mixtures and mixtures of 
alcohols with other polar compoxmds. In addition, systems with limited miscibility have 
been successfully treated with activity coefficient approach. 

2.3 Equation of State Approach 

In the equation of state method, a single equation of state is used to represent all fluid 
phases. From the thermodynamic point of view, this is a potentially more powerful 
approach as it provides a uniform representation of the thermodynamic properties in the 
two-phase and one-phase regions with a single equation. It is not only applicable to wide 
pressure and temperature ranges including critical and supercritical conditions but also 
makes it possible to calculate phase equilibria data. According to the review by 
Kolasinska (1986), depending on the form of equation of state variables, the equation of 
state may be classified into the following types. 

1. Non analytical 

2. Analytical 

a) inspired on a Virial expansion 

b) inspired on the two-term van der Waals form 


Equations of state inspired on the two-term van der Waals form were foimd to be 
most attractive due to their simplicity, low computational costs and reliability. Non- 
analytical types of equations are not presented here. A brief account of the analytical type 
of EOS is presented below. The simplest EOS is the ideal gas equation (PV = nRT) which 
considers that there are no interactions between molecules. 


2.3.1 VirialEOS 


In real gases molecular interactions always exist. This is explained by the Virial 
EOS where PV is expressed as a power series in P (Smith et al. 1996). 

PV = a + bP+cP^ + (2.58) 

PV = a{\ + B'P + C'P^ +....) (2.59) 

Where, a, B\C' are constants at a given temperature for a given species. When a = 
RT, Eq. (2.59) reduces to 

PV 

Z= — = l + B'P+C'P^ + (2.60) 

RT 


or alternatively 


Z 



+ 


C 

V^ 



(2.61) 


Where, Eq. (2.60) or (2.61) is called Virial equation of state, Z is the con^)ressibi]ity 
factor and B, C, D etc., (or B',C',D'etc.,) are called virial coefficients. The virial EOS is 
not used in phase-equilibrium studies as it can describe only vapor phase behavior. 


Equations inspired on the two term van der Waals form 

In general equations based on the van der Waals model are in a form where total 
pressure is the sum of an attractive and a repulsive pressure term (Kolasinska 1986). 

P ~ P (repulsive) T (attractive) (2.62) 

This form is proposed to take into account the attractive and repulsive forces existing 
in the real fluids. Imagine two molecules as hard spheres. When the molecules are at a 
large distance from each other they experience small or no attractive forces. As the 
distance between the molecules is decreased, the attractive forces get larger, until the 



hard spheres come close to each other. Once the molecules come close to each other, 
repulsive forces start dominating until the molecules can get close no further. This 
volume is unavailable for the molecular motion. Van der Waals took the unavailable 
volume into account by introducing a parameter b called the ‘excluded volume’. 


2.3.2 van der Waals Equation of State (vdW) 


When expressed in the form 

P = 


of Eq. (2.62) the van der Waals equation is given by, 

RT a 
v-b 


(2.63) 


repulsive 


RT 

v-b 


(2.64) 



(2.65) 


The values of a and b can be found by regressing vapor pressure data and choosing a 
and b that best fit that data or alternatively critical point conditions are used. At the 
critical point, the first and second derivatives of pressure with respect to volume at 
constant temperature are zero. When these conditions are applied to Eq. (2.63), the values 
of a and b turn out to be 


27R%^ 

64P^ 

b^^ 

8P, 


( 2 . 66 ) 

(2.67) 


2.3.3 Redlich-Kwong Equation of State 

An important modification of the van der Waals equation of state was 
made by Redlich and Kowng (1949), who introduced a temperature dependence and 
slightly different volume dependence in the attractive term. The Redlich-Kwong 
Equation of State (RK EOS) is given by 


RT 


a 


(2.69) 


a = 0.42748 ^ 

Pc 

RT 

b = 0.08664 — £- (2.70) 

Pc 

These modifications improved the accuracy of the equation. As before, the values of a 
and b can be found by fitting vapor pressure data or by using the conditions at the critical 
point. The RK EOS has limited accuracy, and is generally successful only for nearly 
ideal systems. To overcome these deficiencies, many modifications have been proposed 
over the years. Two of these modifications have achieved wide acceptance - the Soave 
and the Peng-Robinson equations of state. 

2.3.4 Soave-Redlich-Kowng Equation of State 


Soave (1972) proposed a modification to the RK Equation incorporating' the 
temperature dependency to the EOS parameter a as 




where 


a(r) = ^ + (0.48 + 1.57 g> - 0.176cu^ 
and CO is the Pitzer (1995) acentric factor defined as 


CO ■■ 


_logl^k±L_i.o 

R. 


(2.71) 

(2.72) 

(2.73) 

(2.74) 

(2.75) 


= (T/T^) = reduced temperature 

Thus, the EOS proposed by Soave is of the form 

J. RT a(r) 
v-b v(v + b) 

When calculating fiigacity coefficients, it is convenient to write the equation in terms of 
the compressibility factor Z. In this form the Soave equation of state is 

-Z^ +z{A-B-B^)-AB = Q (2.76) 

where 

aP 


A = - 


2rp2 


R^T 


(2.77) 


(2.78) 


B 
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bP 

RT 
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(2.79) 


The van der Waals one fluid mixing rules are used for calculating a and b of the mixture. 
Thus the fugacity coefficient of component / in a mixture is given by 
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ln0,. = — (Z-l)-ln(Z-5)+-^ 
b, ’ bRT 
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(2.80) 


2.3.5 Peng-Robinson Equation of State 

In 1976, Peng and Robinson introduced another equation of state, which extended 
some of the ideas presented by Soave. The repulsive part of the van der Waals model 
remained the same, but just as Soave had done, the temperature dependence of the 
attractive term was incorporated into the a term. The denominator of the attractive term 
was also changed to a slightly more complicated expression of volume and the size 
parameter. The Peng-Robinson Equation of State is given by 

RT a(T) 


P = - 


v-b v(v+b) + b(v-b) 

The values of the parameters were obtained by following the method as in SRK EOS, 


(2.81) 


RT 

b = 0.07780 " 

Pc 

(2.82) 

«(rj = 0.45724-^ 

(2.83) 

a(T) = a(T,)a 

(2.84) 

1 1 
=l + K(l-r/) 

(2.85) 

CP 

II 

(2.86) 

i) = 0.37464 + 1 .54226® - 0.26992® ' 

(2.87) 


2.3.6 Peng-Robmson-Stryjek-Vera EOS 

In order to make the Peng-Robinson EOS more widely applicable Stryjck and 
Vera (1986a) modified the temperature and acentric factor dependence of the PR EOS. 
They also extended it to strongly non-ideal mixtures (Stryjek et al, 1986b). They 
introduced a pure component parameter into the expression for k which is slightly 
modified. The PRSV EOS is given by, 

P =— (2.88) 

v-b v(v+b)+b(y-b) 

and the parameters a and b are given by 


RT 

6 = 0.07780 ' 

Pc 

(2.89) 

u(rj = 0.45724^ 

(2.90) 

a{T) = a(T,)a 

(2.91) 

1 1 
=l + ic(l-r,2) 

(2.92) 

where reduced temperature Tr = T/Tc 

In the PRSV equation k term takes the form 

(2.93) 

fc=Xo+K,(l + r;Q(0.7-r,) 

(2.94) 


where jc, is the pure component parameter introduced by Stryjek and Vera and K^is 
given by 

jCo =0.378893 + 1.4897153co -0.17131848®' +0.0196554®' (2.95) 


2.3.7 Quartic Equation Of State 

Kubic (1986) proposed Quartic hard chain EOS and Soave (1990) proposed a new 
Quartic van der Waals type EOS. Shah et al., (1994) introduced the generalized Quartic 
EOS for pure non-polar compounds and extended the same for polar compounds (Shah et 
al, 1996). 

The generalized Quartic EOS is given by. 
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This equation of state is rewritten as quartic in volume as, 
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(2.98) 

(2.99) 

( 2 . 100 ) 
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The hard sphere volumes of the fluids at the critical temperature were fixed as 0.165 
times the critical volume of the fluids. Temperature dependence was incorporated into j3, 
According to Nezbeda and Aim (1984). 


p = p^ exp (-0.03 125 In (r,)- 0.0054 (in (r,))') ' 

(2.102) 

where /3^=0.165F^ 

a = a,a(rj 

(2.103) 

where for < 1 

a(3;) = [l,+ X,(l-V^)+A3(l-77;)%A,(l-V^)'J 

(2.104) 

where for 2^ > 1 


(2.105) 
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II 

(2.106) 

«(r,)=[i+A',(i-r,)]’ 

(2.107) 


where Oc and Cc are the values of a and c at the critical temperature of the fluid, and X 2 . 
X3, X4, Xs. X6 and X7 are the constants determined by regression. The equation of state is 


extended to non-spherical fluids with the introduction of acentric factor cn, as the third 
property to characterize the fluids. The parameters a, c and e, and the constants X 2 
through X? are made quadratic functions of acentric factor. The parameters Oc and Cc, 
related to the acentric factor, are defined as 
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The values of a,o » , c,, 


c^2 > ^ro > » ^r2 Obtained by regression and 


the constants X 2 to X 7 are made functions of co as 

X,=X,^+X,2<^ (2.114) 

for i=2,3,4,5,6,7 


fr 

Extension of the Generalized Quartic EOS to Polar Fluids 


To apply the Quartic equation of sate to polar compounds, the parameters a and c are 
assumed to be functions of the acentric factor and the dipole moment of the fluid. In 
order to express the parameters of the quartic equation in the dimensionless form, the 
reduced dipole moment (jjl) is defined below 

* _ Q.3976p (2.115) 

where n is the dipole moment of the polar fluid in debye units. 

For polar fluids, the regressed coefficients X 2 through X 7 are made linear functions of the 
acentric factor and the quadratic functions of the reduced dipole moment. The regressed 


constants Ur and Cr are made quadratic functions of the acentric factor and the reduced 
dipole moment, and are given by 

X,=Xu-^Xi2+Xi3ix* +Xi4^* (2.116) 
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The generalized Quartic EOS is extended to predict the behavior of the polar flmds. If 
the dipole moment of the fluid is zero, then the entire system is reduced to the non-polar 
fluid system. The values of the regressed coefficients are presented in Table 2. 1 . 


Table 2.1: Quartic Equation of State Constants 


Regressed coefficients 

Regressed coefficients 

Regressed coefficients 

aro 

1.84713 

X 21 

0.14988 

Xs2 

0.57743 

arl 

-0.05218 

X 22 

0.97848 

X 53 

0.41218 

ar2 

1.06446 

X 23 

-0.01390 

X 54 

-0.10676 

ars 

-0.02730 

X 24 

0.02928 

X6, 

0.025SJ 

(lr4 

0.02048 

X}, 

-0.32379 

X62 

-0.02700 

CrO 

1.78336 

X 32 

1.84591 

X63 

0.38327 

Crl 

-1.29690 

X 33 

0.39338 

X64 

-0.09008 

Cr2 

2.78945 

X 34 

-0.25483 

X 71 

-0.77357 

Cr3 

-0.07000 

X 41 

0.14833 

X 72 

-1.45342 

Cr4 

0.01188 

X 42 

-3.46693 

X 73 

-0.04725 

^rO 

0.63189 

X 43 

-0.39170 

X 74 

-0.09669 

eri 

-0.81660 

X 44 

-0.01597 

ko 

1.2865 

^r2 

3.25246 

X 51 

0.11048 

k] 

2.8225 












2.3.8 Advantages and Disadvantages of Equation of state approach 

In the Equation of State approach the definition of fugacity coefficient is used to 
describe both liquid and vapor phases through the application of single equation of state. 
The major advantage of this approach is its applicability over a wide range of pressures 
including critical and supercritical pressures. However, phase equilibrium calculations by 
these methods are rather lengthy except for a simple equation of state. Another problem 
associated with the use of the EOS technique is the insufficiency of information about the 
exact form of the EOS for mixtures, and the inaccuracy of the existing mixing rules for 
mixtures containing polar and hydrogen bonding components. 

2.4 Mixing rules 

Mixing rules used in the prediction of vapor-liquid equilibria of mixtures are as 
important as the equation of state itself Quite often the classical van der Waals one fluid 

mixing rule with one binary interaction parameter was used. However, as the systems 

\ 

become complex in nature, the need for using better mixing rules grows enormously. It is 
highly unlikely that the performance of a new equation of state will be better, unless 
improved mixing rules are used (Tsonopoulas et al., 1985). This led to the proposal of 
many modifications to the existing classical van der Waals type of mixing rules. 

Given below are some of the most widely used mixing rules for a large number 
of equations of state in order to calculate the mixture parameters. 

2.4.1 van der Waals Mixing Rules 

These are also called the classical mixing rules and are given by: 

" = ( 2 . 119 ) 

* j 

^ = (2.120) 

i j 

Where a is the energy parameter and b is the volume parameter, a/,- and bij are the cross 

terms i.e. the interaction terms when (z 5^ J) whereas i mdj represent the components 

(c). 
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Where ky and ly are the binary interaction parameters; ly being zero most of the times. 
When Eq. (2.122) is applied to Eq. (2.120) with ly =0, it reduces to a linear term. 


b = 'Z^ibii 


(2.123) 


The above vdW one fluid mixing rules are used for van der Waals, Redlich-Kwong, 
Soave-Redlich-Kwong, and Peng-Robinson equations of state with a single tnnaty 
interaction parameter ky, (ly = 0). 

2.4.2 Composition Dependent Mixing Rules 


These mixing rules were first proposed by Huron and Vidal (1979) and Stryjek 
and Vera (1986) as modifications to the existing van der Waals mixing rules, to improve 
the vapor-liquid equilibrium calculations. Stryjek and Vera (1986) proposed the 
following mixing rules, 


c c 
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(2.124) 

6 = JxA 

(2.125) 


where the cross terms are same as those for the vdW but, with a change in the binary 
interaction parameter. Two binary interaction parameters ky and kp are used with 
composition dependence for a parameter. Stryjek and Vera (1986) proposed two different 
types of composition dependence for the interaction parameters. 

<^y=^4^^-Xiky-Xjkj) (2.126) 


Eq. (2.126) is a Margules type cross term and 


a.. = Ja-’Ci .. 




I X-k- + X k , 
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Eq. (2. 1 27) is a Van Laar type cross term. . 

The above mentioned mixing rules are used for Peng-Robinson-Stryjek-Vera equation of 
state. 


2.4,3 Conformal Solution van der Waals Mixing Rules 


These mixing rules were proposed by Kwak and Mansoori (1986), and were used 
by Benmekki and Mansoori (1987, 1988). These mixing rules differ depending on the 
EOS used. For the PR EOS, when the temperature dependent term is expanded and re- 
written, 

a{T,(o) = a(TXl + K{\-t‘^)f 


= a(rj(l + 2 k:(1 - + k \\ + T,- 1T°-^)) 

= a{TXK +\f-^a{T;)KX -2jc(x 
= c^-dRT-l4cdRT (2.*128) 

Where 

c = a{f){^+Kf (2.129) 
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b is the same as it is in the PR EOS, leads to the following mixing rules. 
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The cross terms or the combining rules are given by 
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Where, ky, /,y, and rriij are adjustable binary interaction parameters. 

2.4,4 Wong — Sandler Mixing Rule 

Wong and Sandler (1992) developed these mixing rules by equating the Helmholtz free 
energy at infinite pressure obtained by using an equation of state with that obtained from 
a liquid solution model (Activity coefficient model). Orbey and Sandler (1996) used it 
with different EOS to calculate VLE of some highly non-ideal systems. The mixing rules 
are given below: 
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where 


a 


b — ^ is the second virial coefficient obtained by the expansion of the van der 
RT J 


Waals equation in virial form, is the excess Helmholtz; free energy obtained from the 

activity coefficient model, C is a constant specific to the EOS chosen. For the Peng- 
Robinson equation 

C = [ln(V2 -l)J/V2 = -0.62323 
The cross term in Eq. (2.137) is given by 
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Any EOS can use these mixing rules but one has to calculate the excess Helmholtz free 
energy by using any one of the activity coefficient models mentioned previously in this 
chapter. 



CHAPTER 3 


LITERATURE REVIEW 

Many attempts were made in the recent past to determine the VLE of various 
systems and to obtain the model parameters by using experimental VLE data. The work 
done by the Thermodynamic Data Center (TDC) at Warsaw, Poland (Marczynski 1988 ) 
is worth appreciating as they published 6 volumes by collecting the experimental data 
from the world literature. 

The TDC data bank contains experimental data of thermophysical properties for 
pure substances as well as for mixtures. Data for pure substances is primarily collected 
from Thermodynamics Research Center (TRC) at The Texas A&M University, USA, 
while data for binary and ternary mixtures is predominantly from TDC. Vapor-liquid 
equilibrium data is the most essential part of this collection. 

At TDC experimental data is collected on the basis of the open world literature and 
evaluated at the following tests: primary experimental data, selected data, recommended 
data and parameters of selected equation. Selected and recommended data arc prq)ared 
after critical evaluation, which includes correlations and thermodynamic consistency 
tests. Parameters for selected equation are calculated mainly for selected and 
^ recommended data over a given range of temperature. 

The first volume contains selected and correlated data for binary mixtures containing 
Cs - Ce hydrocarbons combined with C5-C20 hydrocarbons. The second volume contains 
the data for binary mixtures containing C7-C14 hydrocarbons combined with CrCu 
hydrocarbons, while the third volume contains selected and correlated data for binary 
mixtures containing CrCi2 alcohols combined with C5-C18 aliphatic hydrocarbons. The 
fourth volume contains the data for binary mixtures containing Ci-Cio alcohols combined 
with C5-G10 non-aliphatic hydrocarbons. The fifth volume deals with the C1-C12 alcohols 
and C3-C10 ethers binary mixtures. The sixth volume contains selected and correlated data 
for binary mixtures containing C3-C18 esters combined with Ci-Ce alcohols and C4-C20 
hydrocarbons. 

With the growing importance of refiigerants, the need for the VLE of the refrigerant 
mixtures was stressed on and many researchers were worked in this area. They also tried 



to obtain the parameters required for various prediction methods by fitting the 
experimental data. Some of those works have been discussed in the subsequent 
paragraphs. Eva Fransson et al. (1992) used a constant -volume cell to determine vapor- 
liquid equilibria by the static method for the binary systems n-pentane + 
chlorodiflouromethane (R22) or + 1,1, -diflouroetane (R152a). They carried out the 
measurements between 303 and 383 K and between 0.3 and 4.4 MPa for different 
compositions of the mixtures. From a knowledge of overall composition, pressure and 
temperature they calculated the vapor and liquid compositions using the Soave-Redlich- 
Kwong equation of state together with an interaction parameter. They concluded that the 
variation of binary interaction parameter with the temperature for these systems is 
negligible. 

The experimental data for vapor liquid equilibrium for the three non-azeotropic 
mixtures HCFC22/CFC114 (CHCIF 2 /C 2 CI 2 F 4 ), HCFC22/HCFC142b (CHCIF 2 /C 2 H 3 CIF 3 ) 
and HCFC22/HFC152a (CHCIF 2 /C 2 H 4 F 2 ) along with the vapor pressures for the ‘pure 
components were determined by Krister et al. (1993). They compared the experimental 
and calculated data and presented the deviations. They used different models for this 
purpose and the models were Lee-kesler and Camahan-Starling-de santis of EOS 
category and activity coefficient model by Wilson. Interaction parameters were fitted 
from experimental data to the mixing rules proposed to the two equations of state and the 
activity coefficient parameters were determined. All three methods studied gave good 
conformity with experimental data for the three mixtures examined. 

Michael Kleiber (1994) obtained the VLE for the binary mixtures propylene-R12, 
propylene-R22, propylene-R152a, propylene-R13, propylene-R23, propylene-R134a, 
propylene-R115, propylene-R13Bl, propylene-R114, propylene-R116, R134a-R12, 
R134a-R152a, R134a-R116, R134a-R12, and R134a-propane at temperatures between 
251 and 298 K and pressure up to 2 MPa. He showed the consistency of the data by a 
maximum-likelihood method. He also correlated the data by generalized forms of various 
EOS such as Peng-Robinson (PR), Peng-Robinson-Stryjeck-vera (PRSV), Redlich- 
kwong-soave (RKS) and Lee-Kesler-Plocker (LKP) along with appropriate mixing rales. 
These EOS were fitted to the experimental data by optimizing the binary interaction 
parameters. Considering the fitting rales. The LKP equation seems to perform worse than 
the others, especially for azeotropic mixtures like propylene-R152a, propylene-R115, 



propylene-R134a, R134a-propane R134a-R12. The PR, RKS and PRSV equations yield 
similar results in each case. The system R134a-R116 was described significantly worse 
than the other mixtures. Only the PRSV equation with the van laar type was able to 
describe it satisfactorily. This mixing rule is, of course, superior to the others, because it 
has two adjustable parameters instead of only one, but usually it gives only slightly better 
results than normal vdW mixing rule. 

Hideo Nishiumi et al. (1995) measured VLB data for the systems- 
chlorodiflouoromethane (HCFC22) - dichlorodifluoromethane (CFC12) over the 
temperature range 348.81 K to 374.14 K and HCFC22 - dichlorotrifluoroethane 
(HCFC123) from 313.48 K to 414.57 K, and vapor pressures and the critical properties of 
their components. Binary interaction parameters of an extended BWR equation of state 
for the systems containing fluorocarbons including the above mentioned systems were 
correlated as quadratic equations of temperature. Using the parameters obtained fix)m the 
VLB measurements they estimated the coefficients of performance (COP) for these 
binary refrigerant mixtures and found that about 17% energy would be saved for the 
system HCFC22-HCFC123 as a mixture-refiigerant in comparison with HCFC22. 

Masanao Kobayashi and Hideo Nishiumi (1998) measured the VLB for the pure, 
binary and ternary systems containing difiuoromethane (HFC32), pentafluoroethane 
(HFC125) and 1,1,1,2-tetrafluoroetane (HFC134a). The binary interaction pamnelers 
' were determined firom the VLB data by using extended BWR equation of state. Using the 
binary parameters, they predicted the VLB data for a ternary system at 333.15 K and 2.43 
MPa, which is in excellent agreement. 

The Isothermal VLB data of the binary system of difiuoromethane (HFC-32) + 1,1,1- 
trifluoroethane (HFC-143a) were obtained in the temperature range from 263.15 K to 
313.15 K by Chang Nyeon Kim and Young Moo Park (2000). They measured the 
temperature, pressure, and compositions of the liquid and vapor phases with a circulation 
type apparatus. The experimental data were correlated with the Camahan-Starling-De 
Santis (CSD), Peng-Robinson, and Redlich-Kwong-Soave equations of state. It is evident 
that the binary interaction parameter for the CSD equation of state increase slightly as the 
temperature increases. On the other hand, the binary interaction parameters for both the 
PR and RKS equations of state decrease as the temperature increases. The binary 
interaction parameter for the RKS equation has the larger temperature dependence than 



those for the CSD and PR equations of state. The average deviation between the 
measured pressures and the calculated results from the CSD equation of state is about 
0.29% and that from the PR equation of state it is about 0.45%. The average deviation 
from the RKS equation of state is about 0.70%. As a result, they concluded that the CSD 
equation of state correlates the experimental data better than the PR and RKS equations 
of state. 

Lim et al. (2000) measured the Isothermal VLB data of the binary mixtures of 
trifluoromethane (HFC-23) + isobutene (R600a) and 1,1,1-trifluoroethane (HFC-143a) + 
isobutene at 283.15 and 293.15 K and at 323.15 and 333.15 K, respectively. They carried 
out the experiment in a circulation-type equilibrium apparatus with the provision for 
measurement of temperature, pressure, and the compositions of the liquid and vapor 
phases. They correlated the experimental data with the Peng-Robinson equation of state 
using the Wong and Sandler mixing rules. When NRTL model was used along with the 
Wong -Sandler mixing rule, the average deviations of pressure are 0.644% for the itFC- 
23 + Isobutane system and 0.30% for the HFC-143a + isobutene system and those of the 
vapor phase compositions are 0.0066 for the HFC-23 + Isobutane system and 0.0058 for 
the HFC-143a + isobutene system. After examining the results they concluded that the 
values calculated with the Peng-Robinson equation of state and Wong sandier mixing 
rules are comparatively in good agreement with experimental data. 



CHAPTER 4 


PREDICTION OF VLE DATA 


The purpose of phase-equilibrium thermodynamics is to predict the conditions (T, P, 
composition) prevailing when two or more phases are in equilibrium. The 
thermodynamic equations which determine the state of equilibrium between viqx>r and 
liquid phases are T' =T''{T is constant), P' =P’'(P is constant), and // = ff .lo find 
the conditions which satisfy these equations, it is necessary to have a method for 
evaluating the fiigacity (/, ) of each component in both vapor and liquid phases. This can 
be done using the basic thermodynamic equation for calculating the fugacity co<^5cient 
(<^/or^;). 


The fugacity coefficient can be obtained by using the relation 
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The solution to the phase-equilibrium problem is provided completely by either one 
of these equations together with an EOS and the equations of phase-equilibrium. Eq. (4.1) 
is used when the EOS is given in volume explicit form (F = Fy(P,T,n^,n 2 ..)) whereas 


Eq. (4.2) is used when the EOS is in the pressure explicit form (P = Py(V,T,n^,n 2 ,..)) . 

Most of the EOS are pressure explicit type. To use an EOS, one must transform the EOS 
in molar form to a total form first. When one is to calculate the fugacity of component i 
in the mixture at a specified T, P, and composition, differentiate the EOS w.r.t rii, and 
substitute in Eq. (4.2) to get the result. 



There are basically four classes of VLE problems: 
BUBLP: Calculate {y,} andP, given {jc,}, T. 
DEWP: Calculate andP, given {yj, T. 
BUBLT: Calculate {y;} and T, given {x,}, P. 
DEWT: Calculate {x,} and T, given {y,}, P. 


The problems of engineering interest generally deal with BUBLP or DEWP. In order 
to correlate the experimental data, choice of the model to calculate the fugacity 
coefficients and the determination of the best or the most representative parameters for 
the model are the most important steps. The models in this chapter are based on the PR 
EOS and Quartic EOS, each combined separately with vdW mixing rules and 
composition dependent two parameter mixing rules. 

4.1 Model 1: PR EOS with van der Waals one fluid mixing rules 

The PR EOS expressed on a total basis, in pressure explicit form, is 

p ^ n,RT «fa(r,6)) 2) 

V-n,b V{V + n,b) + nf{V -n,b) 

where V = vrit and 
fit = total number of moles 


Eq. (4.3) when differentiated on molar basis and substituted in Eq. (4.2) with the use of 
van der Waals one fluid mixing rules (Eqs. (2.119) through (2.123)) gives the following 
expression for the fugacity coefficient. 
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Eq. (4.4) is used to calculate the fugacity coefficient of either liquid or vapor phase, by 
using the appropriate compressibihty factor and the composition. 

4.2 Model H: PR EOS with composition dependent two parameter mixing 
rules 

The PR EOS in the pressure exphcit form expressed on a total basis is, 

p__ n,RT nfa(T,a)) 

V-n,b V{V + n,b) + n,biV-n,b) ^ ^ 

where V = vrit and 

n, = total number of moles 


Eq. (4.3), when used properly in Eq. (4.2), with the use of vdW one fluid and Margules 
type (Stryjek and Vera 1986) mixing rules (Eqs. (2.119), (2.123), and (2.124)) gives the 
following expression for fugacity coefficient. 
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Eq. (4.9) is used to calculate the fugacity coefficients of both vapor and liquid phases. 

4.3 Model III: Quartic EOS with van der Waals one fluid mixing rules 


The Quartic EOS in the pressure explicit form expressed on a total basis is, 

p_ n,RT n^jik^RT n^aV + nXPc 

{V-n,k„p) {V-n,k,py v{V + n,efy-nf,p) ^ ’ 

where V = vn, and 

n, = total number of moles 

Eq. (4.1 1) when differentiated on molar basis and substituted in Eq. (4.2) with the use of 
van der Waals one fluid mixing rules (Eqs. (2.119) through (2.123)) gives the following 
expression for the fugacity coefficient. 
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Eq. (4.12) is used to calculate the fugacity coefficients of both vapor and liquid 
phases. 

4.4 Model IV: Quartic EOS with composition dependent two parameter 
mixing rules 


The Quartic EOS in the pressure explicit form expressed on a total basis is, 

n,RT ^ n^pk,RT n,^aV + n,"kgj3c 
{V-n,kj) {y-n,k,py V{V + n,eyy-n,k,p) 

where V = vn, and 

rit = total number of moles 


Eq. (4.1 1), when used properly in Eq. (4.2), with the use of vdW one fluid and Margules 
type (Stryjek and Vera 1986) mixing rules (Eqs. (2.119), (2.123), and (2.124)) gives the 
following expression for fugacity coefficient. 
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Eq. (4.21) is used to calculate the fiigacity coefficients of both vapor and liquid 
phases. 

4.5 VLE Calculations 

The most important step of calculating the fugacity coefficients (^/ or ) is carried 

out by using the four models described. This requires Z‘ , the liquid phase 
compressibility factor or Z" , the vapor phase compressibility factor. These are obtained 
by solving the EOS which is expressed in terms of the compressibility factor Z. 

For PR EOS, we can obtain the Z'and Z’' by solving the following cubic EOS 
which is expressed in terms of the compressibility factor Z. 


where 
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( 4 . 30 ) 

(iP 


tt 

^ ‘ 

(4.7) 

bP 


B= — 

(4.8) 

RT 


PV 


Z = 

nRT 
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Eq. (4.30) when solved gives three roots. There are two possibilities; either all the three 
roots are real or only one real root (remaining two being complex). When there are three 
real roots, the smallest root is the liquid phase compressibility factor and the largest root 



is the vapor phase compressibiUty factor, intermediate root being meaningless. If it is 
only one root, it is taken either as the liquid phase or vapor phase compressibility factor 
depending on the calculations. 

For Quartic EOS the and Z” are obtained by solving the following quartic EOS 
which is expressed in terms of the compressibility factor Z. 
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On solving Eq. (4.32) one obtains four values of Z, one of which is always negative. 
Among the three positive values of Z, the highest value of Z corresponds to vapor phase 
(Z’^ while the smallest value corresponds to liquid phase (Z^). The intermediate value of 
Z has no physical significance. Substituting Z^ or z' values in Eq. (4.12) or Eq. (4.21) one 


A A 

can get or 0/ 


4.6 Data Reduction 

As mentioned earlier, there are two important steps i) model selection and ii) 
evaluation of the best parameters for the model. This needs reduction of the VLE data. In 
this work, only P-x-y data (given T and xi) is used. Binary interaction parameters in the 
models used are obtained by optimizing an objective function. The objective function 



used in this work is a sum of two functions namely a pressure (OFP) and a vapor phase 
composition (OFY) function (Benmekki et al., 1987). For the optimization of this 
objective function, a function ymins' of the ‘Matlab’, which uses aNelder Mead Simplex 
algorithm is used. 

OF = OFP + OFY (4.37) 

where 
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where, Pexp, Pcai are experimental and calculated pressure respectively, yia^are 


experimental and calculated vapor phase compositions of component 1, and N is the 
number of data points. 

At the optimum values of the binary interaction parameters, the average absolute 
deviations in pressure and composition are calculated fix)m the relations. 


i* 


AP 

P 

AVi 

Jfi 


OFPy 

N ) 

(4.40) 

OFYX^ 

N J 

(4.41) 


The codes employed in this work are given in the Appendix A for all the four models. 

The algorithm used in VLB calculations is given in Fig. (4.1) and Fig. (4.2) and (4.3) 
show the subroutines for PR and Quartic EOS, respectively. 




Fig 4.1. Algorithm for VLB Prediction. 



















Subroutine PR 



Return 


Fig 4.2 Subroutine for Peng-Robinson equation of state. 










Subroutine Quartic 


Input data: 

T,Xi{oryi),P^f,T^^,G)i 


Calculate 
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Calculate the mixture constants depending on the mixing 


Calculate qs, q 2 ,qhqo 
(See Equations 4.33, 4.34, 4.35 and 4.36) 


Solve Tj"' + + ^ 2 ^^ + q^Z 4- — 0 

to obtain 2} or Z' 


Calculate (|./or(|»; using the expression appropriate for the Model 



Fig 4.3 Subroutine for Quartic equation of state 













CHAPTER 5 


RESULTS AND DISCUSSION 


5.1 Results 

The results obtained from the four models are presented in this section. In Ac first 
two models Peng-Robinson EOS was used, and in the next two models Quartic EOS was 
used. The results obtained with these foiu* models are compared to ascertain tfie best 
model. A total of 33 binary systems at different temperatures were studied. These 
systems are classified into 7 different classes, such as alcohol/alcohol, 
alcohol/hydrocarbon, hydrocarbon/hydrocarbon, hydrocarbon / chlorofluorocarbon 
(CFG), hydrocarbon/hydro chlorofluorocarbon (HCFC), hydrocarbon /hydro 

i 

fluorocarbon (HFC) and HFC/HFC. Table (5.1) lists all the systems considered in this 
study. The critical constants, acentric factors and the dipole moments of the compoimds 
(Reid et al. 1987) are listed in Table (5.2). The systems range from non-polar/non-polar 
(hydrocarbon/hydrocarbon) to highly polar/polar (HFC/HFC). Only isothermal VLB 
(BUBLP) data is considered. The objective fimction used for the estimation of the 
optimum binary interaction parameters is a more rigorous one, which is a sum of the 
deviations in both pressure and vapor phase composition. 

The average absolute deviations in pressure and vapor-phase composition along with 
the optimum objective function (OF) of the all systems studied are presented in Table 
(5.3). The binary interaction parameters for all the four models at the optimized objective 
function are presented in Table (5.4). The pressure (P) and vapor-phase composition (y/) 
at the optimum binary interaction parameter values of one system for each class at one 
temperature are presented in Tables (5.5) - (5.11). The minimum percent deviations 
observed for pressure and vapor-composition are 0.12 and 0.35 respectively and the 
corresponding maximum values are 40.92 and 25.26. The results for Methanol/Ethanol 
system are presented in Fig. 5.1. As the predicted results are in close agreement with the 
experimental values the graphical representation for other systems is omitted. 



5.2 Discussion 


In this work thirty three systems of different classes at different temperatures have 
been studied using the PR EOS and Quartic EOS. Two types of mixing rules namely vdW 
mixing rules and composition dependent mixing rule, with Margules type cross term, are 
used. The systems studied in this work range from nonpolar systems (Pentane/Benzene, 
Pentane/Decane, Benzene/hexane, Benzene/Heptane, Benzene/Octane) to highly polar 
systems (R32/R143a). The results obtained indicate that Model 2 i.e. PR EOS with 
composition dependent mixing mles yields better results for most of the systems. For 
some systems, the Model 2 is predicting the VLE data that is as good as the remaining 
models. 

The performances of the models are different for different classes of systems 
considered in this work. For alcohol/alcohol class of systems (Methanol/Propanol, 
Methanol/1 -Butanol, Ethanol/l-butanol), which are polar in nature, the Model 4 i.e. 
Quartic EOS with composition dependent mixing rules predicted the VLE data better 
than the other models. But for Methanol/Ethanol system, the VLE data obtained by the 
Model 4 is as good as that obtained by the Model 2. 

For Alcohol/Hydrocarbon systems, the VLE data obtained by the Model 2 i.e. PR 
EOS with composition dependent mixing rules is better than that by other models. But for 
the systems, which are having 1 -Butanol as of the component (Butanol/Benzene, 
Butanol/Cyclohexane), the VLE data obtained by Model 4 is better than that by other 
models. For hydrocarbon/hydrocarbon systems of class three, which are non polar in 
nature Model 2 is as good as Model 4. But for those systems containing cyclohexane or 
Toluene, which are slightly polar in nature as one of the component of the system 
(cyclohexane/Hexane, Toluene/heptane, Toluene/octane, Toluene/Decane) the VLE data 
obtained by Model 4 is slightly better than Model 2. 

For HC/CFC class systems (Propylene/R12, Propylene/R13, Propylene/RlIS) the 
VLE data obtained by Model 4 is as good as that obtained by Model 2. But in giOMral 
Model 2 is preferable because of its ease in use. For HC/HCFC class of i^nteiiis 
(Propylene/R22, Propylene/R142b) the Model 2 gave better results than other models. 
For HC/HFC class of systems (Propylene/R23, Propylene/R152a) the VLE data obtained 
by the Model 2 is better than that obtained by other models. For HFC/HFC class of 



system (R32/R143a) the Model 2 gave good results. However Model 4 is also equally 
good for this system. 

This study indicates that two parameter composition dependent mixing rules are 
better than the conventional vdW mixing rales for both PR EOS and Quartic EOS. Also 
this study indicates that the well known PR EOS displays better abilities in general as 
compared to the new Quartic EOS for nonpolar systems and refrigerant mixtures. But for 
the polar systems the Quartic EOS gave good results and for polar/nonpolar systems it’s 
results are comparable to the results obtained by PR EOS. 



Table.5.1: List of systems studied and the corresponding temperatures 


S.No 

' SYSTEM 

TEMPERATURE (IQ 

1 

Methanol (1)/ Ethanol (2) 

298.15 

2 

Methanol (1)/ Propanol (2) 

333.17 

3 

Methanol (1)/ 1 -Butanol (2) 

298.15 

4 

Ethanol (1) / 1 -Butanol (2) 

323.15, 343.15, 363.15, 383.15, 403.15 

5 

Methanol (1)/Benzene (2) 

308.15,318.15, 328.15 

6 

Methanol (l)/Toluene (2) 

318.15 

7 

Ethanol (1)/Ben2ene (2) 

298.15,318.15, 323.15,328.22 

8 

Ethanol (l)/Cyclohexane (2) 

273.15, 283.15, 293.15, 298.15, 303.15, 323.15 

9 

Ethanol (l)/Toluene (2) 

308.15,318.15,328.15, 333.15 

10 

1 -Propanol (1)/Benzene (2) 

318.15, 328.15, 348.15 

11 

1 -Propanol (l)/Cyclohexane (2) 

323.01, 328.15, 338.15 

12 

1 -Butanol (1)/Benzene (2) 

298.15,318.15 

13 

1 -Butanol (l)/Cyclohexane (2) 

298.15,318.15, 323.15, 343.15 

14 

1-Butanol (l)/Toluene (2) 

333.31, 343.40, 353.44, 373.15 

15 

Pentane (1)/Benzene (2) 

308.15,313.15,318.15,323.15 

16 

Pentane (1)/Decane (2) 

317.70, 333.70 

17 

Benzene (IZ-Hexane (2) 

298.15, 333.15 

18 

Benzene (l)/Cyclohexane (2) 

283.15, 298.15, 303.15, 313.15, 323.15, 333.15, 

343.15, 403.15,413.15,423.15 

19 

Benzene (1)/Heptane (2) 

293.15, 318.15, 328.15, 333.15, 353.15, 413.15 

20 

Benzene (l)/Toluene (2) 

313.15, 334.15 

21 

Benzene (l)/Octane (2) 

328.15,338.15,348.15 

22 

Cyclohexane (1)/Hexane (2) 

298.15, 308.15, 343.15 

23 

Toluene (1)/Decane (2) 

373.5, 383.6 

24 

Toluene (1)/Heptane (2) 

303.15, 328.15 

25 

Toluene (lyOctane (2) 

333.15 

26 

Propylene (1)/R12 (2) 

258.00, 263.00, 268.00, 273.00, 278.00, 283.00 

27 

Propylene (1)/R13 (2) 

251.00, 273.00 

28 

Propylene (1)/R1 15 (2) 

251.00, 275.00, 298.00 

29 

Propylene (1)-R22 (2) 

258.00, 263.00, 268.00, 273.00, 278.00, 283.00 

30 

Propylene (1)/R142b (2) 

268.00, 278.00, 288.00, 298.00 

31 

Propylene (1)/R23 (2) 

251.00, 265.00 

32 

Pronvlene fl)/R152a (2) 

255.00, 265.00, 275.00, 285.00 





Table 5.1a: Source of Experimental data for the Refrigerant mixtures 


S.No 

System 

Source 

1 

Propylene (1)/R12 (2) 

Michael Kleiber (1994) 

2 

Propylene (1)/R13 (2) 

Michael Kleiber (1994) 

3 

Propylene (1)/R1 15 (2) 

Michael Kleiber (1994) 

4 

Propylene (1)-R22 (2) 

Michael Kleiber (1994) 

5 

Propylene (1)/R142b (2) 

Lim et al. (1999) 

6 

Propylene (1)/R23 (2) 

Michael Kleiber (1994) 

7 

Propylene (1)/R1 52a (2) 

Michael Kleiber (1994) 

8 

R32(l)/R143a(2) 

Chang Nyeon Kim et al. (2000) 




Table 5.2. Properties of the compounds 


H^j 

Compound 

Name 

Tc 

(K) 

IBI 




1 

Methanol 

512.6 

80.9 

118.0 

0.556 

1.70 

2 

Ethanol 

513.9 

61.4 

167.1 

0.644 

1.70 

3 

Propanol 

536.8 

51.7 

219.0 

0.623 

1.70 

4 

Butanol 

563.1 

44.2 

275.0 

0.593 

1.80 

5 

Benzene 

562.2 

48.9 

259.0 

0.212 

0.00 

6 

Toluene 

591.8 

41.0 

316.0 

0.263 

0.40 

7 

Cyclohexane 

553.5 

40.7 

308.0 

0.212 

0.30 

8 

Heptane 

540.3 

27.4 

432.0 

0.349 

0.00 

9 

Hexane 

507.5 

30.1 

370.0 

0.299 

0.00 

10 

Pentane 

469.7 

33.7 

304.0 

0.251 

0.00 

11 

Octane 

568.8 

24.9 

492.0 

0.398 

0.00 

12 

Decane • 

617.7 

21.2 

603.0 

0.489 

0.00 

13 

Propane 

369.8 

42.5 

203.0 

0.153 

0.00 

14 

Isobutane 

408.2 

36.5 

263.0 

0.183 

0.10 

15 

Propylene 

364.9 

46.0 

181.0 

0.144 

0.40 

16 

R23 

299.3 

48.6 

132.7 

0.260 

1.60 

17 

R32 

351.6 

58.3 

120.8 

0.271 

2.00 

18 

R143a 

346.3 

37.6 

194.0 

0.251 

2.30 

19 

R152a 

386.7 

45.0 

181.0 

0.256 

2.30 

20 

R12 

385.0 

41.4 

216.7 

0.204 

0.54 

21 

R22 

369.3 

49.7 

165.6 

0.221 

1.40 

22 

R13 

302.0 

38.7 

180.4 

0.198 

0.50 

23 

R142b 

409.6 

43.3 

231.0 

0.251 

2.10 

24 

R115 
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0.34 
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Table 5,3 Optimized average absolute deviations of vapor phase composition and pressure and Objective Function (OF) 
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0.0054 
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0.0067 

0.0076 

0.0056 

0.0027 
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0.0197 

0.0210 

0.0214 

0.0225 

0.0376 

0.0345 

0.0309 

0.0307 

0.0241 

0.0299 

0.0226 

0.0273 

0.0216 

0.0211 

0.0180 

0.0208 

0.0201 

0.0233 

0.0181 

0.0154 
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0.0461 

1 

0.0185 

! 0.0188 

0.0204 

0.0195 

0.0375 

0.0346 

0.0312 

0.0306 

0.0393 

0.0373 

0.0316 

0.0229 

0.0277 

0.0221 

0.0211 

0.0417 

0.0406 

0.0465 

0.0395 

0.0455 

AY/Y 

0.0169 

0.0211 

0.0181 

0.0228 

0.0178 

0.0072 

0.0108 

0.0132 

0.0183 

0.0527 

0.0413 

0.0146 

0.0141 

0.0115 

0.0125 

0.0513 

0.0526 

0.0615 

0.0563 

0.0526 

PRM2 

OF 

0.0003 

0.0007 

0.0003 

0.0010 

0.0017 

0.0023 

0.0067 

0.0120 

0.0014 

0.0081 

0.0035 

0.0034 

0.0016 

0.0026 

0.0016 

0.0023 

0.0032 

0.0076 

0.0047 

0.0026 

§ 

0.0045 

0.0072 

0.0055 

0.0072 

0.0139 

0.0140 

0.0196 

0.0295 

0.0054 

0.0133 
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|g 

268.00 j 

273.00 

278.00 

283.00 

268.00 

278.00 

288.00 

298.00 

251.00 

265.00 

255.00 

265.00 

275.00 

285.00 

263.15 

273.15 

283.15 

293.15 

303.15 

313.15 

SYSTEM 

Propylene(l )-Rl 42 b (2) 

Propylene ( l )- R 23(2) 

Propylene ( l)-Rl 52 a (2) 

R 32( l )- R 143 a (2) 
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Table 5.4. Optinmm binary interaction parameters of the systems 
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Table 5.4. Optimum binary interaction parameters of the systems (continued.) 
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Table 5.4. Optimum binary interaction parameters of the systems (continued.) 
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Table 5.4. Optimum binary interaction parameters of the systems (continued.) 
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Table 5.4. Optimum.binary interaction parameters of the systems (continued.) 
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Table 5.5: VLE Data for Methanol (1) - Ethanol (2) 
Temperature: 298.15 K 
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Table 5.5: VLE Data for Methanol (1) — Ethanol (2) (continued.) 

Temperature: 298.15 K 
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Table 5.6: VLE Bata for Ethanol (1) - Benzene (2) 
Temperature: 298.15 K 
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Table 5.6i VLE Data for Ethanol (1) — Benzene (2) (continued.) 

Temperature: 298.15 K 


Model 4 

P (Kpa) 

13.1082 
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16.3878 
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7.3988 

0.0328 

1 

00 
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0.0000 
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0.3861 
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0.1077 

o 

Model 3 

P(Kpa) 

13.1082 

14.4220 

15.6423 
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16.1197 

16.2928 

16.3208 

16.3111 
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ZZZVO 

Experimental Data 

P (Kpa) 

12.7700 
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16.9300 
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16.4800 
16.0900 
15.2400 

14.7300 
14.4400 
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7.9700 


Ayjy 

Objective Function 
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Table 5.7: VLE Data for Benzene (1) - Toluene (2) 
Temperature: 334,15 K 
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Table 5.7: VLE Data for Benzene (1) - Toluene (2) (continued.) 

Temperature: 334.15 K 


Model 4 

P(Kpa) 

19.3405 

21.5416 

23.3581 

25.2382 

27.5006 

29.8035 

31.0574 

34.2186 

35.9803 

39.7162 

42.1819 

46.9688 

48.7728 

51.0887 

53.0305 

55.1188 

ZLOOO 


! 

o 

o 

o 

o 


0.0000 

0.1539 

0.2599 

0.3543 

0.4516 

0.5362 

0.5772 

0.6678 

0.7117 

0.7918 

0.8366 

0.9094 

0.9327 

0.9598 

0.9801 

1.0000 


0.0032 

1 Experimental Data 1 Model 3 

j 1 

q. 

Cd 

Ph 

j 19.3405 
! 21.5297 
23.3492 
25.2434 
27.5358 
29.8807 
31.1609 
34.3935 
36.1946 
39.9999 
42.4925 
47.2646 
49.0347 
51.2816 
53.1420 
55.1188 

0.0078 


91000 

0.0000 

0.1531 

0.2587 

0.3527 

0.4495 

0.5335 

0.5741 

0.6637 

0.7069 

0.7861 
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0.9041 

0.9281 

0.9562 

0.9781 

1.0000 


£900 0 

19.2900 

21.5720 

23.3290 

25.1690 

27.5210 

29.9360 

31.2100 

34.3260 

36.2040 

39.9920 

42.3280 

46.9170 

48.5570 

50.6450 

52.3490 

54.0590 

d/dv 

$ 

Objective Function 


0.0000 

0.1535 

0.2589 

0.3516 

0.4519 

0.5390 

0.5802 

0.6683 

0.7128 

0.7949 

0.8372 
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0.9332 

0.9574 

0.9810 

1.0000 


1 1 

i f 

0.0000 

0.0622 
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0.3003 

0.3368 

0.4288 

0.4799 

0.5873 

0.6571 

0.7891 

0.8374 

0.8981 

0.9478 

1.0000 





Table 5.8: VLB Data for Propylene (1) - R12 (2) 
Temperature: 258.00 K 









Table 5.8: VLE Data for Propylene (1) - R12 (2) (continued.) 

Temperature: 258.00 K 


<D 

O 


m 

o 


Ph 


cd 

Ph 

w 

Oh 




oJ 

P. 

Nd 

Ph 




169.5218 

228.7176 

253.6015 

278.9108 

298.3044 

322.6322 

345.4399 

359.4898 

376.4215 

0.0305 


o 

K 

00 

to 

00 

CO 

00 

o 

o 



o 

O 

:3: 


o 

o 

o 

o 

o 


o\ 

o 

T* 

CO 

xt* 

CM 

T" 

CD 


o 


o 

o 


to 

CO 

N- 

00 

00 

CD 

o 


o 

o 

d 

d 

d 

d 

d 

d 

d 

T~ 


d 

00 

to 

to 

00 

CD 

K 



to 


■1 

T*" 

00 

CM 

to 


00 

00 

o 

X— 



CM 

(0) 

00 

00 

to 

to 

00 

CM 

CM 

0\ 


ID 

xf 

xf 

q 

to 

o 

00 

00 

xf 

o 

CO 

■ 

o> 

00 

o6 

00 

CO 

o6 

to 

CD 

CD 

d 

■ 

CD 

CM 

to 

K 

CD 

CM 

xj- 

lO 





CM 

CM 

CM 

CM 

00 

00 

CO 

CO 


■ 

o 

K 

K 

00 

to 

CD 

CD 


o 


00 

o 

O) 

00 

00 

CD 

00 

CD 

CD 

o 


r-* 

o 

O 

00 

xf 

T- 

o 

00 

CO 

o 


o 

o 


to 

CD 

N- 

00 

00 

q 

q 


o 

d 

d 

d 

d 

d 

d 

d 

d 

T" 


d 

o 

o 

o 

o 

o 

o 

o 

o 

o 



o 

o 

o 

o 

o 

o 

o 

o 

o 



o 

o 

o 

o 

o 

o 

o 

o 

o 



00 

00 

h- 

q 

CD 

o 

q 

o 




’T— 

to 

00 

d 


d 


CD 

CM 



00 

00 

to 

00 

CD 

CM 

CO 

Xt 

CD 



■r" 

CM 

CM 

CM 

CM 

oo 

oo 

CO 

CO 



o 

T- 

to 

h- 

K 

oo 

CD 

CD 

o 

Cl. 


o 

00 

O) 

to 

to 

xf 

o 

o 

o 



o 

o 

CM 

00 


T~ 

CD 

Xf 

o 

< 

q 


to 

q 

b- 

00 

00 

q 

q 



d 

d 

d 

d 

d 

d 

d 

d 




o 

T** 

'•ct 

xf 

to 

h- 

CO 

CD 

o 



o 

0) 

CM 

LO 

xt 

oo 

00 

o 

o 



o 

00 

to 

N- 

N- 

o 

CM 

CD 




q 

CM 

00 


q 

N. 

cq 


q 



d 

d 

d 

d 

d 

d 

d 

Vm/ 





o^ 

o 


o 

o 


Ph 

o 

> 

‘p 

o 

<D 

o 


VO 






Table 5.9: VLE Data for Propylene (1) - R22 (2) 
Temperature; 258.00 K 
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Table 5.9: VLE Data for Propylene (1) - R22 (2) (continued.) 

Temperature: 258.00 K 




LO 

CO 

CM 

CD 

00 

CO 

to 

to 

to 


■ 



cd 

<J> 

CM 

CO 

CM 

hi 

CO 

00 

00 

CD 

CM 

o 

On 

■ 




00 


CD 



CD 

00 


▼-H 




o 

CD 

CM 


d 

cd 

to 

cd 

o 

1 



Ph 

C3) 

C\J 

CM 

CO 

'SJ- 

CO 

CO 

CO 

N. 

CO 

K 

CO 

N. 

CO 

N- 

CO 

d 

■ 

xt 

X) 











■ 

O 

Cr 















O 

to 

CO 

ID 

oo 

xl- 

CO 

o 

■ 

CM 

tn 

d 



O 

CM 

s 

CD 

CO 

N. 

to 

o 





O 

IT" 

Tj- 

CD 

CO 

CM 

o 


H 




o 

CM 


CO 

CO 

h- 

00 

o 

■ 

o 




CD 

CD 

d 

d 

d 

d 

d 


1 

d 




to 

T- 

CD 

CD 


CD 

CO 

q 


to 





Cd 

CD 

o 

N. 

CD 

N. 

to 

T“ 

1 





T** 

CD 

K 

CO 

CO 

to 

CM 





& 

00 

CD 

CO 

CD 

00 

CM 

Tt- 

r-H 



m 

cd 

CM 

CD 

to 

00 

CM 

xt 

d 

q 



Ph 

CD 

CM 

xt 

CO 

CO 

r*^ 

CO 

N. 


d 




CM 

CO 

CO 

CO 

CO 

CO 

CO 



r-H 














XJ 

o 












o 












q 



o 



CO 

CD 

to 

T- 

o 



o 



o 


K 

CO 

o 

T— 

CD 

o 


00 




o 

T— 

to 

to 

o 

N- 

CM 

o 


r-H 



q 

CM 


CO 

K 

h- 

00 

q 


o 




o 

d 

d 

d 

d 

d 

d 



d 




o 

o 

o 

o 

o 

o 

o 

o 






o 

o 

o 

o 

o 

o 

o 

o 





Cu 

o 

o 

o 

o 

o 

o 

o 

o 




B 


q 

CO 

T- 

N. 

cd 

cd 

q 

q 

cd 

q 

Tt 

CM 




Ph 

CD 

CM 


CO 

CO 

CO 

CO 

CO 



q 


CM 

CO 

CO 

CO 

CO 

CO 

CO 

CO 



o 

P 












• rH 

o 

-j 












a 

B 

Cj 


O 

N. 

o 

CO 

T— 

CM 

00 

o 



s 

CD 


O 

CO 

CO 

00 

CO 

CO 

CO 

o 



<D 


o 

T-“ 

xjf 

CO 

00 

to 


o 

< 

> 

B 

q 

q 


q 

q 

q 

cq 

q 


• rH 

4-H 

G 

CD 


CD 

d 

d 

d 

d 

d 

d 

T- 



o 

CD 


























O 

PQ 


o 

h- 

o 

o 

CO 

to 

CO 

o 






o 

T" 

CO 

K 

CO 

CO 

CO 

o 






o 

iX) 

00 

o 

CO 


T“ 

o 





K ! 

o 

T" 

CO 

q 

CO 

h: 

q 

q 






CD 

d 

d 

d 

d 

d 

d 

T" 













Table 5.10: VLEBata for Propylene (1) - R152a (2) 

Temperature: 255.00 K 
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Table 5.11: VLE Data for R32 (1) - R143a (2) 
Temperature; 263.15 K 
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Table 5.11: VLE Data for R32 (1) - R143a (2) (continued) 

Temperature: 263.15 K 
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Figure 5.1 P-x-y diagram for Methanol (l)/Ethanol (2) system at 298.15 k 


CHAPTER 6 


CONCLUSIONS 

Based on the present study it is evident that the VLB data predicted by the PR EOS 
with composition dependent mixing rule of Margules type with two binary interaction 
parameters is in good agreement with the experimentally obtained VLB data. This 
conclusion is valid for most of the systems, especially the nonpolar systems and 
refrigerant mixtures. 

Tlic Quartic EOS is expected to predict the VLB data better in view of the presence 
of dipole moment parameter and its complexity when compared to PR EOS. But this 
appears to be true only for highly non polar mixtures only. From the present study it can 
be concluded that the two binary interaction parameter mixing mles predicts the VLB 
data much better tlum vdW single binary interaction parameter mixing rule. This is true 
for both PR EOS and Quartic EOS. 

b'urther studies can be carried out for the prediction of VLB data of a large 
number of systems using both these EOS to judge which one is superior. As a 
continuation of the work presented, it is also suggested that two parameter mixing mles 
vi/,. the Van Eaar type, be used. 
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Appendix A 


A1 Program Code for VLE data by Model I 

% Matlab commands contained in the file “Modell.m”. 

function ModellQ; %Main flmction i.e. calls other functions to calculate optimum 

%bmary interaction parameters and VLE data corresponding to 
%that values. 

min=jBnin('pr',-0. 1 ,0.23) %”min” is the vector containing the final interaction parameters 

% where the objective function becomes minimum. ”pr” is the 
%function for VLE.”finin” subfunction minimizing a function 
%one variable using golden section search algorithm and the two 
%values are the boundaries in which one tries to find the 
% minima. 

% Control is passed onto the file “pr.m”. “pr.m” contains the following commands. 

Function OF=pr(kij) %OF is the objective function value which is passed to the 

%“finin” function and “kij” is the value passed by “finin”. 

global I; 

1=2; %Initialization of the required data. 

j=l ;STATEL=0;STATEV=1 ; 

n=8; %Benzene&Heptane ’ 

Xexpt=[0.0000, 0.3516, 0.5252,0.5633, 0.6083, 0.6513, 0.6731, 1.0000]; % Vector with experimental 

Yexpt=[0.0000,0.4748,0.6479,0.6671,0.7192,0.7351,0.7536,1.0000]; %Vector with experimental 
' Pexpt=[298.5,384.2, 41 8.6,425.0,430.9,436.0,440.6,473.2]; %experimentarF 

T=413.15; % temperature 

Tc=[562.2, 540.3]; %Critical temperature of components 1 and 2. 

Pc=[48.9,27.4]; %Critical pressure of components 1 and 2. 

vc=[0.259,0.432]; %Critical volumes of components 1 and 2. 

w=[0.2 12,0.349]; % Acentric factors of components 1 and 2 

for j=l:n % Loop changing experimental ‘x’ values 

x(l)=Xexpt(j); 
x(2)=(l-x(l)); 

P=0; 

psat(2)=Pexpt(l); 

psat(l)=Pexpt(n); 

P=x(l)*psat(l)+x(2)*psat(2); %Initial guess of total pressure 

sumyl=0; 

P1=0; 


while(sumyl~=l) 



if(abs(P-P 1)<1 e-05) 
break; 
end 

fori=l:I 

k(i)=psat(i)/P; 

y(i)=k(i)*x(i); %Initial guess for vapor composition ‘y 

end 
t=i; 

[phi]= pil (T,P,x,Tc,Pc,w,STATEL,kij); %Control passes to “pi 1 .m” where fugacity coefficient 

%for liquid phase is calculated. The code for “pil.m” is 
% given at the bottom. 


phil=phi; 

condl=l;cond2=l; 

wliile(l) 

if((condl=0)&(cond2=0)) 

break; 

end 

t=l; 

[phi]=pil(T,P,y,Tc,Pc,w,STATEV,kij); 

phiv=phi; 
sumyl=0; 
for i=l :I 

y 1 (i)=(phil(i)*x(i))/phiv(i); 
sumyl=sumyl+yl(i); 
end 

for i=l :I 
if(i==l) 

if(abs(y(i)-yl(i))<1.0e-06) 
condl=0; 

end 
end 

if(i=2) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
cond2=0; 
end 
end 

y(i)=yl(i); 
end 
end 


% Fugacity coefficient for the vapor phase is 
% calculated using “pil.m” 


% y is calculated fi-om fugacity coefficients. 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 1. 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 2. 

%Change the value of ‘y’ to new. 

% end of inner “while” loop 


P1=P; 

P=P*sumyl; 

end 

YcalO’)=y(l); 

Pthe(j)=P; 

end 


%New guess value for pressure 
% end of outer “while” loop 
%Put calculated values of ‘y’ into vector. 
%Put calculated values of pressure into vector. 
% end of outer “for”loop 


OF=0;OFY=0;OFP=0; 



if((Ycal(l)=0)&(Yexpt(l)=0)) 

OFY=0; 

else 

0FY=0FY+(((Yexpt(l)-Ycal(l))/Yexpt{l))'^2); 

end 

for i=2:n 

OFY=OFY+(((Yexpt(i)-Ycal(i))A?'expt(i))^2); %Objective function of %composition. 
end 

for i=l :n 

OFP=OFP+(((Pexpt(i)-Pthe(i))/Pexpt(i))'^2); %Objective function of pressure, 
end 

OF=OFY+OFP %Overall objective function. 

%Code for function “pil.m” is given below. 


function [phi]=pil(T,P,x,Tc,Pc,w, state, ki) %calculates fugacity coefficient. 
R = 8.314; 


for i=l :2 

s(i)=0.37464+(l .54226*w(i))-(0.26992*(w(i)^2)); % 

alpha(i)=((l+s(i)*(l-sqrt(T/Tc(i))))^2); % PR EOS parameters. 

ai(i)=0.45724*(((R*Tc(i))^2)*alpha(i))/(Pc(i)* 1 00000); % 
bi(i)=(0.07780*R*Tc(i))/(Pc(i)* 100000); % 


end 

a=amix(x,ai,ki); 

b=bmix(x,bi); 

A=(a*P* 1 000)/((R*T)^2); 

B=(b*P*1000)/(R*T); 

p=(B-l); 

q=(A-(3.0*(B^2))-(2.0*B)); 
n=((B^3)+(B^2)-(A*B)); 
ro=[l p q r]; 

Z=Toots(ro); 


% EOS parameter ‘a’ for mixture 
%EOS parameter ‘6’ for mixture 


%-(l-B) 

% coefficients of the cubic equation of compressibility factor. 
% 

% 

% Roots (Compressibility factor) are found out. 


Big=0;Small=0; 

if(state=l) 

if((imag(Z(l))==0)&(real(Z(l))>=0)) 

Big=real(Z(l)); 

end 

if((imag(Z(2))=0)&(real(Z(2))>=0)) 

if((real(Z(2)))>(real(Z(l)))) 

Big=real(Z(2)); 

end 

end 

if((imag(Z(3))==0)«&(real(Z(3))>=0)) 

if(Big<(real(Z(3)))); 

Big=(real(Z(3))); 

end 

end 

Zl=Big; 


% 

% 

% 

%AiTange the compressibility factors 
%as vapor and liquid phase 
% compressibility factors 
% 

% 

% 

% 

% 

% 

% 

% 

% 



end 


% 
% 

if(state==0) % 

if((imag(Z(l))=0)&(real(Z(l))>=0)) % 

Small=real(Z(l)); % 

End % 

if((imag(Z(2))=0)&(real(Z(2))>=0)) % 

if((real(Z(2)))<(real(Z(l)))) % 

Small=real(Z(2)); % 

end ' % 

end % 

if((imag(Z(3))=0)&(real(Z(3))>=0)) % 

if(Small>(real(Z(3)))) % 

Small=(real(Z(3))); % 

elseif(Small=0) % 

Small=real(Z(3)); % 

end % 

end % 

Zl=Small; % 

end % 


al 2=(x( 1 )* ai( 1 )+(x(2)*(sqrt(ai( 1 )* ai(2)))*( 1 .0-ki))); 

a21=(x(2)*ai(2)+(x(l)*(sqrt(ai(l)*ai(2)))*(1.0-ki))); 

phi(l)=exp(((bi(l)/b)*(Zl-l))-(log(Zl-B))+(a/(2*sqrt(2)*b*R*T))*((bi(l)^)- %Fugacity 

((2*al2/a)))*(log((Zl+B*(l+sqrt(2.0)))/(Zl+B*(l-sqrt(2.0)))))); % coefficients 

% of both the 

phi(2)=exp(((bi(2)/b)*(Zl-l))-(log(Zl-B))+(a/(2*sqrt(2)*b*R*T))*((bi(2)^)- % components. 
((2*a21/a)))*(log((Zl+B*(l+sqrt(2.0)))/(Zl+B*(l-sqrt(2.0)))))); % 

% Code for calculating mixture parameters ‘a’ and ‘b’ “amix.m” and “bmix.m” 

function a = amix(x,ai,ki) 

a =((x(ir2)*ai(l)+(x(2)^2)*ai(2j+2.0*x(l)*x(2)*sqrt(ai(l)*ai(2))*(l .0-ki)); 

function b=bmix(x,bi) 
b==x( 1 )*bi(l )+x(2)*bi(2); 


A2 Program Code for VLB data by Model II 

% Matlab commands contained in the file “Model2.m”. 

function Model20; %Main function i.e. calls other flmctions to calculate optimum 
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%binary interaction parameters and VLE data corresponding to 
%that values. 

kij= [0,0] 

min=f mins ( ' prsv ' , ki j ) %”min” is the vector containing the final interaction parameters 

% where the objective function becomes minimum. ”prsv” is the 
%function for VLE.”fmin” subfunction minimizing a function 
%one variable using golden section search algorithm and the two 
%values are the boundaries in which one tries to find the 
% minima. 

%Control is passed onto the file “prsv.m”. “prsv.m” contains the following commands. 

function OF=prsv (ki j ) %OF is the objective function value which is passed to the 

%“jBmin” function and “kij” is the value passed by “finin”. 

global I; 

1=2; %Initialization of the required data. 

j=l ;STATEL=0;STATEV=1 ; 
n=8; %Benzene&Heptane 

Xexpt=[0.0000, 0.3516,0.5252, 0.5633,0.6083, 0.6513, 0.6731,1-0000]; % Vector with experimental 'x 
Yexpt=[0.0000,0.4748, 0.6479,0.6671,0.7192,0.7351,0.7536,1.0000]; %Vector with experimental 'y 
Pexpt=[298.5,384.2,418.6,425.0, 430.9,436.0,440.6,473.2]; %experimentarP' 

T=4 13.15; % temperature 

Tc=[562.2,540.3]; %Critical temperature of components 1 and 2. 

Pc=[48.9,27.4]; %Critical pressure of components 1 and 2. 

vc=[0.259,0.432]; %Critical volumes of components 1 and 2. 

w=[0.212,0.349]; % Acentric factors of components 1 and 2 

for j=l :n % Loop changing experimental ‘x’ values 

x(l)=Xexpt(j); 
x(2)=(l-x(l)); 

P=0; 

psat(2)=Pexpt(l); 

psat(l)=Pexpt(n); 

P=x(l )*psat(l)+x(2)*psat(2); %Initial guess of total pressure 

sumyl=0; 

P1=0; 

while(sumyl~=l) 

if(abs(P-Pl)<le-05) 

break; 

end 

for i=l :I 
k(i)=psat(i)/P; 

y(i)=k(i)*x(i); %Initial guess for vapor composition ‘y’. 

end 

t=l; 


83 



[phi]= prsvpil(T,P,x,Tc,Pc,w,STATEL,kij);%Control passes to “prsvpiLm” where fugacity 

%coefficient for liquid phase is calculated. The code fc 
%“prsvpil.m”is given at the bottom. 

phil=phi; 

condl=l;cond2=l; 

wlhle(l) 

if((condl~0)&(cond2=0)) 

break; 

end 


t=l; 

[phi]=prsvpil(T,P,y,Tc,Pc,w,STATEV,kij); % Fugacity coefficient for the vapor phase is 

% calculated using “prsvpil.m” 

phiv=phi; 
sumyl=0; 
for i=l :I 

yl(i)=(phil(i)*x(i))/phiv(i); % y' is calculated from fugacity coefficients, 

sumy l=sumy 1+y 1 (i); 
end 

for i=l :I 
if(i=l) 

if(abs(y(i)-y 1 (i))<l .Oe-06) % Check if old and new values for ‘y’ arc within 

condl=0; , % tolerable limit for component 1. 

end 
end 


if(i==2) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
cond2=0; 
end 
end 

y(i)=yl(i); 

end 

end 


P1=P; 

P=P*sumyl; 

d 

YcalG)=y(l); 

Pthe(j)=P; 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 2. 


%Change the value of ‘y’ to new. 

% end of inner “while” loop 

%New guess value for pressme 
% end of outer “while” loop 
%Put calculated values of ‘y’ into vector. 
%Put calculated values of pressure into vector. 
% end of outer “for’loop 


OF=0;OFY=0;OFP=0; 

if((Ycal(l)==0)&(Yexpt(l)=0)) 

OFY=0; 

else 

OFY=OFY+(((Yexpt(l)-Ycal(l))A?'expt(l))^2); 


end 


for i=2:n 

OFY=OFY+(((Yexpt(i)-Ycal(i))/Yexpt(i))^2); %Objective function of %composition. 
end 
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fori=l:n 

0FP=0FP+(((Pexpt(i)-Pthe(i))/Pexpt(i))^2); %Objective function of pressure, 
end 

OF=OFY+OFP 


%Overall objective function. 


%Code for function “prsvpil.m” is given below. 


function [phi]=pil(T,P,x,Tc,Pc,w, state, ki) %calculates fligacity coefficient. 

R = 8.314; 
for i=l :2 

s(i)=0.37464+(1.54226*w(i))-(0.26992*(w(i)^2)); % 

alpha(i)=((l+s(i)*(l-sqrt(T/Tc(i))))'^2); % PR EOS parameters. 

ai(i)=0.45724*(((R*Tc(i))^2)*alpha(i))/(Pc(i)* 100000); % 
bi(i)=(0.07780*R*Tc(i))/(Pc(i)*100000); % 

end 

% EOS parameter ‘a’ for mixture 
%EOS parameter ‘b’ for mixture 


a=amix(x,ai,ki); 
b=bmix(x,bi); 

A=(a*P* 1 000)/((R*T)^2); 
B=(b*P*1000)/(R*T); 

P=(B-1); 

q=(A-(3.0*(B^2))-(2.0*B)); 
r=((B^3)+(B^2)-(A*B)); 
ro=[l p q r]; 

Z=roots(ro); 

Big=0;Small=0; 

if(state=l) 

if((imag(Z(l))==0)&(real(Z(l))>=0)) 

Big=real(Z(l)); 

end 

if((imag(Z(2))=0)&(real(Z(2))>=0)) 

if((real(Z(2)))>(real(Z(l)))) 

' Big=real(Z(2)); 

end 

end 

if((imag(Z(3))=0)&(real(Z(3))>=0)) 

if(Big<(real(Z(3)))); 

Big=(real(Z(3))); 

end 

end 

Zl=Big; 

end 

if(state==0) 

if((imag(Z(l))==0)&(real(Z(l))>=0)) 

Small=real(Z(l)); 

End 

if((imag(Z(2))=0)&(real(Z(2))>=0)) 

if((real(Z(2)))<(real(Z(l)))) 


%-(l-B) 

% coefficients of the cubic equation of compressibility factor. 
% 

% 

% Roots (Compressibility factor) are found out. 


% 

% 

% 

%AiTange the compressibility factors 
%as vapor and liquid phase 
% compressibility factors 
% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 
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Small=Teal(Z(2)); 

end 0 /^ 

end 0 /^ 

if((imag(Z(3))==0)&(real(Z(3))>=0)) % 

if(Small>(real(Z(3)))) o/^ 

Small=(real(Z(3))); o/^ 

elseif(Small— 0) % 

Small=real(Z(3)); o/^ 

end % 

end % 

Zl=SmaU; 

end % 


al2=(x(l)*ai(l)+x(2)*sqrt(ai(l)*ai(2))*(1.0-ki(l)*x(l)- 

ki( 2 )*x( 2 ))+x(l)*(x( 2 )^ 2 )*sqrt(ai(l)*ai( 2 ))*(ld( 2 )-ki(l))); 

a21=(x(2)*ai(2)+x(l)*sqrt(ai(l)*ai(2))*(1.0-ki(l)=^x(l)- 

ki(2)*x(2))-(x(2)*(x(l)^2)*sqrt(ai(l)*ai(2))*(ki(l)-ki(2))); 

phi( 1 )=exp(((bi( 1 )/b)*(Z 1 - 1 ))-(Iog(Zl -B))+(a/(2*sqrt(2)*b*R*T))*((bi(l)^)- 
(2*^al2/a))*(log((Zl+B*(l+sqrt(2.0)))/(Zl+B*(l-sqrt(2.0)))))); 

phi(2) cxp(((bi(2)/b)*(Zl-^))-(log(Zl-B))+(a/(2*sqrt(2)♦b*R*T))*((bi(2)^)- 
(2*a21/a))*(log((ZH-B*(H-sqrt(2.0)))/(Zl+B*(l-sqrt(2.0)))))); 

% Code for calculating mixture parameters ‘a’ and ‘b’ “prsvamix.m” and “prsvbmix.m” 

function a » pravamix(x,ai,ki) 

a * ( ( (x (1) ■'2) *ai (1) ) + ( (x(2) '‘2) *ai (2)) + (2.0*x(l) *x{2) *sgrt (ai (1) *ai (2) ) * 

(1.0- ki (1) *X (1) -ki (2) *x(2) ) ) ) ; 


funct: ion bsprsvbmix (x,bi) 
b-=x(l) *b.i (1) +x{2) *bi (2) ; 


A3 Program Code for VLB data by Model III 

% Matlab commands contained in the file “Modelql.m”. 

fimclion Modclq I (); %Main function i.e. calls other functions to calculate optimum 

%binaiy interaction parameters and VLB data corresponding to 
%that values. 

min rmin('qpr’,-(). 1 ,0.23) %”min” is the vector containing the final interaction parameters 

% where the objective fimction becomes minimum. ”qpr” is the 
%function for VLE.”finin” subfimction minimizing a function 
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%one variable using golden section search algorithm and the two 
%values are the boundaries in which one tries to find the 
% minima. 

%Control is passed onto the file “qpr.m”. “qpr.m” contains the following commands. 

Function OF=qpr(kij) %OF is the objective fimction value which is passed to the 

%“finin” function and “kij” is the value passed by “finin”. 

global I; 

1=2; %Initialization of the required data. 

j=l ;STATEL=0;STATEV=1 ; 
n=8; %Benzene&Heptane 

Xexpt=[0.0000, 0.3516,0.5252, 0.5633, 0.6083,0.6513, 0.6731, 1.0000]; % Vector with experimental 'x 
Yexpt=[0.0000,0.4748,0.6479,0.6671,0.7192,0.735 1,0.7536,1 .0000]; %Vector with experimental 'y 
Pexpt=[298.5,384.2,41 8.6,425.0,430.9,436.0,440.6,473.2]; %experimentarP' 

T=413.15; 

Tc=[562.2,540.3]; 

Pc=[48.9,27.4]; 
vc=[0.259,0.432]; 
mu=[0.0,0.0]; 

for j=l:n 
x(l)=Xexpt(j); 
x(2)=(l-x(l)); 

P=0; 

psat(2)=Pexpt(l); 
psat(l)=Pexpt(n); 

P=x(l)*psat(l)+x(2)*psat(2); %Initial guess of total pressure 

surayl=0; 

P1=0; 

' while(sumyl~=l) 

if(abs(P-Pl)<le-05) 

break; 

end 

fori=l:I 

k(i)=psat(i)/P; 

y(i)=k(i)*x(i); %Initial guess for vapor composition ‘y’. 

end 

t=l; 

[phi]= qpil(TJ>,x,Tc,Pc,mu,STATEL,kij);%Control passes to “qpil.m” where fugaaty coefficien 

%for liquid phase is calculated. The code for “qpi 1 .m” 
% given at the bottom. 


% temperature 

%Critical temperature of components 1 and 2. 

%Critical pressure of components 1 and 2. 

%Critical volumes of components 1 and 2. 

%Dipole moment of components 1 and 2 

% Loop changing experimental ‘x’ values 


phil=phi; 

condl=l ;cond2=l ; 
while(l) 

if((cond 1 =0)&(cond2=0)) 
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break; 

end 


t=l; 

[phi]=qpil(T,P,y,Tc,Pc,mu,STATEV,kij); % Fugacity coefficient for the vapor phase is 

% calculated using “qpil.m” 

phiv=phi; 
sumyl=0; 
for i=l:I 

yl(i)=(phil(i)*x(i))/phiv(i); % y' is calculated from fugacity coefficients, 

sumy l=sumy 1+y 1 (i); 
end 


fori=l:I 


if(i=l) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
condl=0; 
end 
end 

if(i=2) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
cond2=0; 
end 
end 

y(i)=yl(i); 

end 

end 


P1=P; 


P=P*sumyl; 

end 

Ycal(j)=y(l); 

PtheO')=P; 

end 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 1 . 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 2. 


%Change the value of ‘y’ to new. 

% end of inner “while” loop 

%New guess value for pressure 

% end of outer “while” loop 

%Put calculated values of ‘y’ into vector. 

%Put calculated values of pressure into vector. 
% end of outer “for”loop 


OF=0;OFY=0;OFP=0; 

if((Ycal(l)=0)&(Yexpt(l)=0)) 

OFY=0; 

else 


OFY=OFY+(((Yexpt(l)-Ycal(l))A"expt(l))'^2); 

end 


for i=2:n 

OFY=OFY+(((Yexpt(i)-Ycal(i))A^expt(i))''2); %Objective function of %composition. 
end 


for i=l :n 

OFP=OFP+(((Pexpt(i)-Pthe(i))/Pexpt(i))^2); %Objective function of pressure, 
end 

OF=OFY+OFP %Overall objective function. 


% Code for function ’'qpil.m" is given below. 
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function [phi]=qpil(T,P,x,Tc, mu, vc, state, ki) %calculates fiigacity coefficient. 

R= 8.314/1000; 
k0=1.2865; 
kl=2.8225; 
phi=[0,0]; 
for i=l :2 

mus(i)=(0.3976*mu(i)/(sqrt(R*Tc(i)*vc(i)))); 

X2(i)=0.14988+0.97848*w(i)-0.01390*mus(i)+0.02928*mus(i)*mus(i); 

X3(i)=-0.32379+1.84591*w(i)+0.39338*mus(i)-0.25483*mus(i)*mus(i); 

X4(i)=0.14833-3.46693*w(i)-0.39170*mus(i)-0.01597*mus(i)*mus(i); 

X7(i)=-0.77357-1.45342*w(i)-0.04725*mus(i)-0.09669*mus(i)*mus(i); 

bei(i)=0.165*(vc(i))*((exp(-0.03125*(log(T/Tc(i)))-0.0054*((log(T/Tc(i)))^2)))'"3); 

ei(i)=(0.63189)*(l-0.81660*w(i)+3.25246*w(i)*w(i))*(vc(i)); 

ai(i)=((l+X2(i)*(l-sqrt(T/Tc(i)))+X3(i)*((l-sqrt(T/Tc(i)))"'2)+X4(i)* 

((l-sqrt(T/Tc(i)))^3))'^2)*(1.84713)*(l-0.05218*w(i)+l-06446*w(i)*w(i) 

-0.02730*mus(i)+0.02048*mus(i)*mus(i))*(R*vc(i)*Tc(i)); 
ci(i)=((l+X7(i)*(l-sqrt(T/Tc(i))))"^2)*(l .78336)*(1-1 .29690*w(i)+ 2.78945 *w(i)Mi) 
+0.07000*mus(i)+0.01188*mus(i)*mus(i))*(R*vc(i)*Tc(i)); 
end 

a=qamix(x,ai,ki); % EOS parameter 'a' for mixture 

be=qbemix(x,bei); %EOS parameter 'beta' for mixture; 

c=qcmix(x,ci,ki); %EOS parameter 'c' for mixture 

e=qemix(x,ei); %EOS parameter 'e' for mixture 

q3=(P/(R?'T))*(-2*kO*be+e-(R*T/P)); % Coefficients of the 

q2=((P/(R*T))^2)*(((R*T/P)*(be*(kO-kl)-e))+kO*be*(kO*be-2*e)+a/P); % quartic equation 

ql=((P/(R*T))^3)*(e*(((kO*be)^2)+(R*T*be/P)*(kO-kl))+(kO*be*(c-a)/P)); % of compressibility 
qO=((P/(R*T))^4)*(-c*((kO*be)^2)/P); % factor 

ro=[l q3 q2 ql qO]; 

Z=roots(ro) ; % Roots (Compressibility factor) are found out. 

Big=0;Small=0; 

if(state=l) 

'if((imag(Z(l))==0)«&(real(Z(l))>=0)) 

Big=real(Z(l)); 

end %Arrange the compressibility factors 

if((imag(Z(2))==0)&(real(Z(2))>=0)) %as vapor and liquid phase 

if((real(Z(2)))>(real(Z(l)))) % compressibility factors 

Big=real(Z(2)); 
end 
end 

if((imag(Z(3))==0)&(real(Z(3))>=0)) 

if(Big<(real(Z(3)))) 

Big=(real(Z(3))); 

end 

end 

if((imag(Z(4))=0)&(real(Z(4))>=0)) 

if(Big<(real(Z(4)))) 

Big=(real(Z(4))); 
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end 

end 


Zl=Big; 

end 


if(state=0) 

if((imag(Z(l))=0)&(real(Z(l))>=0)) 

Small=Teal(Z(l)); 

end 

if((imag(Z(2))=0)&(real(Z(2))>=0)) 

if((real(Z(2)))<(real(Z(l)))) 

Small=real(Z(2)); 

end 

end 

if((imag(Z(3))=0)&(real(Z(3))>=0)) 

if(Small>(real(Z(3)))) 

Small=(real(Z(3))) ; 
end 
end 

if((iniag(Z(4))=0)&(real(Z(4))>=0)) 

if(Small>(real(Z(4)))) 

Small=(real(Z(4))); 

elseif(Small=0) 

Small=real(Z(4)); 

end 

end 

Zl=Small; 

end 

vl=Zl*R*T/P; 

El=(P*e/(R*T)); 

E2=((P*kO*be)/(R*T)); 
aij(l)=x(l)*ai(l)+x(2)*sqrt(ai(l)*ai(2))*(l-ki); 
aij(2)=x(2)*ai(2)+x(l)*sqrt(ai(l)*ai(2))*(l-ki); 
cij ( 1 )=x( 1 )*ci( 1 )+x(2)*sqrt(ci( 1 )*ci(2))* ( 1 -ki); 
cij (2)=x(2)*ci(2)+x( 1 )*sqrt(ci(l )*ci(2))*(l -ki); 


fori=l;2 

LI (i)=((kO*bei(i))/(vl -(kO*be))); 

L2(i)=((kl*bei(i)*vl)/((vl-(k0*be))'^2)); 

L3(i)=((kl*be)/(vl-(k0*be))); 

L4(i)=-((l/(R*T))*(a/(kO*be+e))*((kO*bei(i)/(vl-(kO*be)))+(ei(i)/(vl+e)))); 

L5(i)=((l/(R*T))*(c/(kO*be+e))*(((kO*bei(i)*ei(i))/(e*(vl+e)))-((kO*bei(i)/(vl-kO*be))))); 

L(i)=Ll(i)+L2(i)+L3(i)+L4(i)+L5(i); 

M 1 (i)=(2*aij (i)/ ((kO*be)+e)); 

M2(i)=-((a*(k0*bei(i)+ei(i))+c*(ei(i)-(k0*be*k0*bei(i)/e)))/((k0*be+e)^2)); 

M3(i)=((kO*be*2*cij(i)+kO*bei(i)*c)/(kO*be*(kO*be+e))); 
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M4(i)=-(k0*bei(i)*c/(k0*be*e)); 

M(i)=Ml(i)+M2(i)H-M3(i)+M4(i); 

Nl(i)=-((2*aij(i))/((k0*be)+e)); 

N2(i)=((a*(k0*bei(i)+ei(i))+c*(ei(i)-(k0*be*k0*bei(i)/e)))/((k0*be+e)'^2)); 

N3(i)=((k0*be*2*cij(i)+k0*bei(i)*c)/(e*(k0*be+e))); 

N4(i)=-(ei(i)*c/(e*e)); 

N(i)=Nl(i)+N2(i)+N3(i)+N4(i); 

Ql(i)=-((k0*be*2*cij(i)+k0*bei(i)*c)/(k0*be*e)); 

Q2(i)=(c*(bei(i)*e+be*ei(i))/(be*e*e)); 

Q(i)=Ql(i)+Q2(i); 

phi(i)=exp((L(i)+(l/(R*T))*(Q(i)*(log(Zl))+N(i)*(log(Zl+El))+M(i)*(log(Zl-E2)))-(log(Zl-E2)))) 

end 

% Code for calculating mixture parameters ‘a%‘beta’,’c’ and ‘e’ 
“qamix.m”,“qbemix.m”,”qcmix.m” and “qemix.m” 
function a = qamix(x,ai,ki) 

a=((x(ir2)*ai(l)+(x(2)^2)*ai(2)+2.0*x(l)*x(2)*sqrt(ai(l)*ai(2))*(1.0-ki)); 

function be=qbemix(x,bei) 
be=x(l)*bei(l)+x(2)*bei(2); 

function c = qcmix(x,ci,ki) 

c=((x(ir2)*ci(l)+(x(2)"^2)*ci(2)+2.0*x(l)*x(2)*sqrt(ci(l)*ci(2))*(1.0-ki)); 

function e = qemix(x,ei) 
e=x(l)*ei(l)+x(2)*ei(2); 


A4 Program Code for VLB data by Model IV 

% Matlab commands contained in the file ‘‘Modelq2.m”. 

%Main function i.e. calls other functions to calculate optimum 
%binary interaction parameters and VLE data corresponding to 
%that values. 

%”min” is the vector containing the final interaction parameters 
% where the objective function becomes minimum. ”qsv” is the 
%function for VLE.”finin” subfunction minimizing a function 
%one variable using golden section search algorithm and the two 
%values are the boundaries in which one tries to find the 
% minima. 

%Control is passed onto the file “qsv.m”. “qsv.m” contains the following commands. 

Function OF=qpr(kij) %OF is the objective function value which is passed to the 


function Modelq20; 
kij=[0,0] 

min=fniins(’qsv',kij) 
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%“fi3iin” function and “kij” is the value passed by “finin”. 

global I; 

1=2; %Initialization of the required data. 

j=l ;STATEL=0;STATEV=1 ; 
n=8; %Benzene&Heptane 

Xexpt=[0.0000, 0.3516, 0.5252, 0.5633, 0.6083, 0.6513, 0.6731,1.0000]; % Vector with experimental 
Yexpt=[0.0000,0.4748,0.6479,0.6671,0.7192,0.7351,0.7536,1.0000]; %Vector with experimental 
Pexpt=[298.5,384.2,418.6,425.0,430.9,436.0,440.6,473.2]; YoexperimentalT' 


T=413.15; 

Tc=[562.2,540.3]; 

Pc=[48.9,27.4]; 

vc=[0.259,0.432]; 

mu=[0.0,0.0]; 


% temperature 

%Critical temperature of components 1 and 2. 
%Critical pressure of components 1 and 2. 
%Critical volumes of components 1 and 2. 
%Dipole moment of components 1 and 2 


for j=l :n % Loop changing experimental ‘x’ values 

x(l)=XexptO); 
x(2)=(l-x(l)); 

P=0; 

psat(2)=Pexpt(l); 

psat(l)=Pexpt(n); 

P=x(l)*psat(l)+x(2)*psat(2); %Initial guess of total pressure 

sumyl=0; 

P1=0; 


while(sumyl~=l) 

if(abs(P-Pl)<le-05) 

break; 

end 

fori=l:I 

k(i)=psat(i)/P; 

y(i)=k(i)*x(i); %Initial guess for vapor composition ‘y’. 

> end 
t=l * 

[phi]= qpil(T,P,x,Tc,Pc,mu,STATEL,kij); %Control passes to “qsvpil.m” where fugacity 

% coefficient for liquid phase is calculated. The code 
% for “qsvpil.m” is given at the bottom. 

phil=piii; 

condl=l;cond2=l; 

while(l) 

if((cond 1 =0)&(cond2=0)) 
break; 
end 
t=l; 

[phi]=qpil(T,P,y,Tc,Pc,mu,STATEV,kij); % Fugacity coefficient for the vapor phase is 

% calculated using “qsvpil .m” 

phiv=phi; 
sumyl=0; 
for i=l:I 



y 1 (i)=(phil(i)*x(i))/phiv(i); 
sumyl=sumyl+yl(i); 
end 

fori=l:I 

if(i==l) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
condl=0; 
end 
end 


% y' is calculated from fugacity coefficients. 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 1. 


if(i=2) 

if(abs(y(i)-y 1 (i))<l .Oe-06) 
con^=0; 
end 
end 

y(i)=yl(i); 

end 

end 


P1=P; 


P=P*sumyl; 

end 

YcalO)=y(l); 

PtheG)=P; 

end 


OF=0;OFY=0;OFP=0; 


% Check if old and new values for ‘y’ are within 
% tolerable limit for component 2. 


%Change the value of ‘y’ to new. 

% end of inner ‘ Vhile” loop 

%New guess value for pressure 
% end of outer “while” loop 
%Put calculated values of ‘y’ into vector. 
%Put calculated values of pressure into vector. 
% end of outer “for”loop 


if((Ycal(l)=0)&(Yexpt(l)==0)) 

OFY=0; 

else 


OFY=OFY+(((Yexpt(l)-Ycal(l))A"expt(l))^2); 

end 

for i=2:n 

OFY=OFY+(((Yexpt(i)-Ycal(i))A?'expt(i))^2); %Objective function of %composition. 
end 


for i=l :n 

OFP=OFP+(((Pexpt(i)-Pthe(i))/Pexpt(i))^2); %Obj ective fiinction of pressure, 
end 

OF=OFYh-OFP %Overall objective function. 


% Code for function "qsvpil.m" is given below. 

function [phi]=qpil(T,P,x,Tc, mu, vc, state, ki) %calculates fugacity coefficient. 

R = 8.314/1000; 

k0=1.2865; 

kl=2.8225; 

phi=[0,0]; 

for i=l :2 

mus(i)=(0.3976*mu(i)/(sqrt(R*Tc(i)*vc(i)))); 

X2(i)=0.14988+0.97848*w(i)-0.01390*mus(i)+0.02928*mus(i)*mus(i); 
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A3(l)=-u.32i/y+1.845yi'^w(l)+u.3y338’‘mus(l>u.254«j5‘^mus(l)‘^mus(l); 

X4(i)=0.14833-3.46693*w(i)-0.39170*mus(i)-0.01597*mus(i)*mus(i); 

X7(i)~0.77357-1.45342*w(i)-0.04725*mus(i)-0.09669*mus(i)*mus(i); 

bei(i)=0. 1 65*(vc(i))*((exp(-0.03 125*(log(T/Tc(i)))-0.0054*((log(T/Tc(i)))^2)))^3); 

ei(i)=(0.63 1 89)*(l-0.8 1 660*w(i)+3.25246*w(i)*w(i))*(vc(i)); 

ai(i)=((l+X2(i)*(l-sqrt(T/Tc(i)))+X3(i)*((l-sqrt(T/Tc(i)))^2)+X4(i)* 

((l-sqrt(T/Tc(i)))"^3))'^2)*(1.84713)*(l-0.05218*w(i)+l-06446*w(i)*w(i) 

-0.02730*mus(i)+0.02048*mus(i)*mus(i))*(R*vc(i)*Tc(i)); 
ci(i)=((l+X7(i)*(l-sqrt(T/Tc(i))))^2)*(l .78336)*(1-1 .29690*w(i)+ 2.78945*w(i)*w(i) 
+0.07000*mus(i)+0.01 1 88*mus(i)*mus(i))*(R*vc(i)*Tc(i)); 
end 

a=qamix(x,ai,ki); % EOS parameter 'a' for mixture 

be=qbemix(x,bei); %EOS pai ameter 'beta' for mixture; 

c=qcmix(x,ci,ki); %EOS parameter 'c' for mixture 

e=qemix(x,ei); %EOS parameter 'e' for mixture 

q3=(P/(R*T))*(-2*kO*be+e-(R*T/P)); % Coefficients of the 

q2=((P/(R*T))'^2)*(((R*T/P)*(be*(kO-kl)-e))+kO*be*(kO*be-2*e)+a/P); % quartic equation 

ql=((P/(R*T)r3)*(e*(((kO*be)^2)+(R*T*be/P)*(kO-kl))+(kO*be*(c-a)/P));% of compressibility 
qO=((P/(R*T))M)*(-c*((kO*be)^2)/P); % factor 

ro=[l q3 q2 ql qO]; 

Z=roots(ro) ; % Roots (Compressibility factor) are found out. 

Big=0;Small=0; 

if(state=l) 

if((imag(Z(l))=0)&(real(Z(l))>=0)) 

Big=real(Z(l)); 

end %Arrange the compressibility factors 

if((imag(Z(2))==0)&(real(Z(2))>=0)) %as vapor and liquid phase 

if((real(Z(2)))>(real(Z(l)))) % compressibility factors 

Big=real(Z(2)); 
end 
end 

if((imag(Z(3))=0)&(real(Z(3))>=0)) 

' if(Big<(real(Z(3)))) 

Big=(real(Z(3))); 

end 

end 

if((imag(Z(4))=0)&(real(Z(4))>=0)) 

if(Big<(real(Z(4)))) 

Big=(real(Z(4))); 

end 

end 

Zl=Big; 

end 

if(state=0) 

if((imag(Z( 1 ))==0)&(real(Z( 1 ))>=0)) 

Small=Teal(Z(l)); 
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if((imag(Z(2))=0)&(real(Z(2))>=0)) 

if((real(Z(2)))<(real(Z(l)))) 

Small=real(Z(2)); 

end 

end 

if((imag(Z(3))=0)&(real(Z(3))>=0)) 

if(Small>(real(Z(3)))) 

Small=(real(Z(3))); 

end 

end 

if((imag(Z(4))=0)&(real(Z(4))>=0)) 

if(Small>(real(Z(4)))) 

Small=(real(Z(4))); 

elseif(Small=0) 

Sniall=real(Z(4)); 

end 

end 

Zl=Small; 

end 

vl=Zl*R*T/P; 

El=(P*e/(R*T)); 

E2=((P*kO*be)/(R*T)); 

aij(l)=x(l)*ai(l)+x(2)*sqrt(ai(l)*ai(2))*(l-x(l)*ki(l)- 

x(2)*ki(2))+x(l)*(x(2)^2)*sqrt(ai(l)*ai(2))*(ki(l)-ki(2)); 

aij(2)=x(2)*ai(2)+x(l)*sqrt(ai(l)*ai(2))*(l-x(l)*ki(l)- 

x(2)*ki(2))+(x(l)^2)*x(2)*sqrt(ai(l)*ai(2))*(ki(2)-ki(l)); 

cij(l)=x(l)*ci(l)+x(2)*sqrt(ci(l)*ci(2))*(l-x(l)*ki(l)- 

x(2)*ki(2))+x(l)*(x(2)^2)*sqrt(ci(l)*ci(2))*(ki(l)-ki(2)); 

cij (2)=x(2)*ci(2)+x( 1 )*sqrt(ci( 1 ) *ci(2))*( 1 -x( 1 )*ki( 1 )- 
x(2)''ki(2))+(x(l)^2)*x(2)*sqrt(ci(l)*ci(2))*(ki(2)-ki(l)); 

for i=l :2 

LI (i)=((kO*bei(i))/(vl -(kO*be))); 

L2(i)=(0cl*bei(i)*vl)/((vl-(k0*be))'^2)); 

L3(i)=((kl *be)/(vl -(kO*be))); 

L4(i)=-((l/(R*T))*(a/(kO*be+e))*((kO*bei(i)/(vl-(kO*be)))+(ei(i)/(vl+e)))); 

L5(i)=((l/(R*T))*(c/(kO*be+e))*(((kO*bei(i)*ei(i))/(e*(vl+e)))-((kO*bei(i)/(vl-kO*be))))); 

L(i)=Ll(i)+L2(i)+L3(i)+L4(i)+L5(i); 

M 1 (i)=(2*aij(i)/((k0*be)+e)); 

M2(i)=-((a*(k0*bei(i)+ei(i))+c*(ei(i)-(k0*be*k0*bei(i)/e)))/((k0*be+e)^2)); 

M3(i)=((k0*be*2*cij(i)+k0*bei(i)*c)/(k0*be*(k0*be+e))); 

M4(i)=-(k0*bei(i)*c/(k0*be*e)); 

M(i)=Ml (i)+M2(i)+M3(i)+M4(i); 
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N1 (i)=-((2*ay(i))/((k0*be)+e)); 

N2(i)=((a*(k0*bei(i)+ei(i))+c*(ei(i)-(k0*be*k0*bei(i)/e)))/((k0*be+e)^2)); 

N3(i)=((k0*be*2*cij(i)+k0*bei(i)*c)/(e*(k0*be+e))); 

N4(i)=-(ei(i)*c/(e*e)); 

N(i)=Nl (i)+N2(i)+N3(i)+N4(i); 

Ql(i)— ((k0*be*2*cij(i)+k0*bei(i)*c)/(k0*be*e)); 

Q2(i)=(c*(bei(i)*e+be*ei(i))/(be*e*e)); 

Q(i)=Ql(i)+Q2(i); 

phi(i)=exp((L(i)+(l/(R*T))*(Q(i)*(log(Zl))+N(i)*(log(Zl+El))+M(i)*(log(Zl-E2)))-(log(Zl-E2)))) 

end 


% Code for calculating mixture parameters ‘a%‘beta’,’c’ and ‘e’ 
“qsvamix.m”,“qsvbemix.m”,”qsvcmix.m” and “qsvemix.m” 

function a = qsvamix(x,ai,ki) 

a=((x(ir2)*ai(l)+(x(2)^2)*ai(2)+2.0*x(l)*x(2)*sqrt(ai(l)*ai(2))*(1.0-x(l)*ki(l)-x(2)*ki(2))); 


flinction be=qsvbemix(x,bei) 
be=x(l)*bei(l)+x(2)*bei(2); 

function c = qsvcmix(x,ci,ki) 

c=((x(ir2)*ci(l)+(x(2r2)*ci(2)+2.0*x(l)*x(2)*sqrt(ci(l)*ci(2))*(1.0-x(l)*ki(l)-x(2)*ki(2))); 

I* 

flinction e = qsvemix(x,ei) 
e=x(l )*ei(l )+x(2)*ei(2); 
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